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ABSTRACT 

A  second-order  epsilon  method  is  developed  for  trajec- 
tory optimization  problems.   The  method  is  applied  to  sever- 
al aircraft  and  missile  performance  and  air  combat  maneuver- 
ing problems.   Heavy  emphasis  is  placed  on  the  realistic 
modeling  of  the  flight  vehicle's  motion  and  maneuvering 
limitations . 

The  proposed  optimization  technique,  which  is  an  exten- 
sion of  Balakrishnan' s  epsilon  method,  uses  either  the  full 
second-order  Newton-Raphson  method  or  the  "modified"  Newton- 
Raphson  method  to  minimize  the  epsilon  functional.   The  full 
Newton-Raphson  method  exhibits  terminal  convergence  charac- 
teristics superior  to  the  "modified"  method,  whereas  the 
"modified"  method  is  generally  superior  in  the  initial 
stages  of  a  problem.   An  algorithm  is  developed  which  uses 
both  techniques  in  a  complementary  way. 

A  new  penalty  functional  which  has  desirable  theoretical 
properties  and  exhibits  excellent  computational  behavior  is 
introduced  to  treat  state  and  control  inequality  constraints. 
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I.   INTRODUCTION 

The  objective  of  the  research  reported  on  herein  is  to 
develop  a  method  of  solving  realistic  problems  in  aircraft 
and  missile  performance  optimization.   Optimization  problems 
of  this  type  have  been  the  subject  of  considerable  research 
[Refs.  1,  2,  3,  H ,  5>  6,  and  7].   The  mathematical  models 
used  in  these  references  are  the  products  of  many  simplifi- 
cations and  assumptions.   Typically,  the  degree  of  simplifi- 
cation used  to  render  these  problems  solvable  by  some  opti- 
mization technique  is  such  that  the  solutions  obtained  are 
of  limited  practical  value.   This  is  particularly  true  in 
the  modeling  of  aircraft  maneuvering  limitations,  such  as 
aerodynamic  stall,  maximum  structural  load  factor,  and 
placard  Mach  number,  which  require  the  use  of  multiple 
state  and  control  inequality  constraints.   Since  these 
limitations  play  an  extremely  important  role  in  maneuvering 
flight  and  air  combat,  they  must  be  modeled  as  accurately 
as  possible.   It  is,  therefore,  imperative  that  the  optimi- 
zation technique  used  to  solve  the  problems  posed  herein 
be  capable  of  handling  state  and  control  inequality 
constraints  with  relative  ease. 

Balakrishnan' s  epsilon  method  is  an  attractive  optimiza- 
tion technique  because  of  the  natural  manner  in  which  state 
and  control  inequality  constraints  are  introduced.   The 
epsilon  method  is  a  penalty  method  in  which  terms  are  added 


to  the  performance  measure  to  penalize  deviations  from  the 
state  equations  written  as  equality  constraints.   Likewise, 
state  and  control  inequality  constraints  may  be  treated  by 
the  addition  of  appropriate  penalty  terms  to  the  performance 
measure.   The  resulting  augmented  performance  measure  is 
minimized  by  an  appropriate  algorithm  for  solving  unconstrained 
optimization  problems. 

The  optimization  technique  used  most  successfully  in 
the  literature  with  the  epsilon  method  is  a  "modified" 
Newton-Raphson  technique,  hereafter  referred  to  as  the  MNR 
technique.   In  this  method  certain  second-order  terms  pre- 
sent in  the  full  Newton-Raphson  formulation  are  neglected. 
It  is  argued  [Ref.  8]  that  since  computer  storage  and  time 
requirements  to  compute  these  terms  are  large,  and  since 
satisfactory  results  can  be  obtained  over  a  large  class  of 
problems  without  the  terms,  their  inclusion  is  not  justified. 
For  these  reasons  the  full  Newton-Raphson  formulation,  here- 
after referred  to  as  the  FNR  technique,  has  not  been 
previously  utilized  with  the  epsilon  method. 

Difficulties  were  experienced  by  the  author,  however, 
in  applying  the  MNR  technique  to  problems  of  the  type  formu- 
lated in  this  dissertation.   The  MNR  technique  was  not 
effective  in  problems  with  state  equations  and  multiple 
inequality  constraints  resulting  from  a  realistic  modeling 
of  the  flight  vehicle's  motion  and  maneuvering  limitations. 

For  this  reason  the  FNR  method  was  investigated  and 
found  to  be  feasible  in  terms  of  computational  storage  and 


time  requirements.   The  PNR  method  exhibits  terminal  con- 
vergence characteristics  superior  to  the  MNR  method  although 
the  MNR  method  is  generally  superior  in  starting  a  problem. 
Problems  not  solvable  by  the  MNR  method  alone  were  solved 
by  an  algorithm  which  uses  both  methods  in  a  complementary 
way.   The  power  of  the  FNR  technique  close  to  the  minimum 
can  also  be  used  to  advantage  to  obtain  a  family  of  optimal 
trajectories  for  different  end  conditions.   The  optimal 
trajectory  for  one  set  of  end  conditions  is  used  as  a  first 
guess  for  the  optimal  trajectory  for  a  neighboring  set  of 
end  conditions. 

Several  simplified  problems  in  aircraft  performance 
optimization  were  attempted  initially  to  gain  experience 
with  the  epsilon  method.   In  these  problems  inequality 
constraints  were  treated  by  using  interior  penalty  func- 
tionals  of  the  type  recommended  by  Fiacco,  McCormick,  and 
Jones  [Refs.  9  and  10].   Computational  results  were  unsatis- 
factory.  Difficulties  were  experienced  in  keeping  the  con- 
strained state  or  control  completely  admissible;  a  require- 
ment for  the  success  of  an  algorithm  with  this  type  of 
penalty  functional.   To  alleviate  this  difficulty  a  new 
penalty  functional  for  inequality  constraints  is  introduced 
which  exhibits  excellent  computational  behavior.   The  pro- 
posed functional  has  performed  well  in  computation  with  up 
to  eight  inequality  constraints  represented  in  a  single 
problem. 
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The  thesis  is  divided  into  eight  sections.   In  Section 
II  the  epsilon  method  is  presented.   The  FNR  and  MNR  tech- 
niques are  derived  and  discussed.   The  effectiveness  of 
the  FNR  method  as  opposed  to  the  MNR  method  is  demonstrated 
by  a  scalar  example.   Finally,  the  computational  experience 
gained  with  both  methods  in  solving  realistic  performance 
problems  is  presented.   In  Section  III  the  method  of  treating 
state  and  control  inequality  constraints  is  presented.   The 
author's  experience  with  interior  penalty  methods  is  related 
and  a  new  penalty  functional  is  proposed.   Computational 
experience  with  the  new  penalty  functional  is  related  and, 
finally,  several  desirable  theoretical  properties  of  the 
new  penalty  functional  are  presented.   In  Section  IV  the 
algorithm  developed  for  minimizing  the  epsilon  functional 
by  either  the  MNR  or  FNR  methods  is  presented.   In  Sections 
V,  VI,  and  VII  three  aircraft  and  missile  performance  opti- 
mization problems  are  solved.   These  problems  are  pertinent 
and  realistic  in  their  operational  applicability.   The 
three-degree-of-freedom  models  are  the  same  as  those  used 
in  basic  aircraft  performance  analysis.   Finally,  the 
summary  and  conclusions  are  presented  in  Section  VIII. 


II.   THE  EPSILON  METHOD 

This  section  describes  the  epsilon  method  and  reviews 
the  significant  contributions  of  other  investigators.   The 
full  Newton-Raphson  (FNR)  equations  for  minimizing  the  aug- 
mented performance  measure  are  derived  and  compared  to  the 
modified  Newton-Raphson  (MNR)  equations  published  elsewhere 
[Refs.  8,  11,  and  12]. 

A.   DESCRIPTION  OF  THE  EPSILON  METHOD 

1.  Statement  of  the  Problem 

A  dynamic  system  characterized  by  the  nonlinear 
state  equations 

x(t)  =  f[x(t),u(t),t]  (2.1) 

is  to  be  controlled  to  minimize  the  performance  measure 

J(x,u)  =  h[x(T),T]  +   /  g[x(t),u(t),t]dt     (2.2) 

where  x(t)  is  an  n  x  1  state  vector  and  u(t)  is  an  I   x  1 
control  vector.   State  and  control  inequality  constraints 
are  omitted  for  the  present.   In  Section  III  the  inclusion 
of  these  constraints  is  discussed  in  detail. 

2.  The  Augmented  Performance  Measure 

In  the  epsilon  method  as  proposed  by  Balakrishnan 
[Refs.  13  and  14],  the  performance  measure  (2.2)  Is  augmented 
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by  a  penalty  functional  which  involves  a  weighted  integral 
of  the  Euclidean  norm  of  the  state  equations  written  as 
equality  constraints.   The  augmented  performance  measure  is 


T  2 

J  (x,u,e)  =  J(x,u)  +  i  f   ||x  -  f(x,u,t)||  dt    (2 


3) 


=  J(x,u)  +  i  J  (x,u).  (2.4) 

#v    *w  t-      O    -^    <w 


The  weighting  factor  e  is  a  positive  quantity. 

3.   Behavior  as  e  -»  0 

As  e  is  reduced,  the  penalty  term  J   is  more  heavily 

s 

weighted,  thereby  placing  greater  emphasis  on  satisfying 
the  state  equations.   Balakrishnan  [Refs.  13  and  14]  and 
Taylor  [Ref.  11]  have  shown  that  under  appropriate  assump- 
tions as  e  ■*■   0,  the  epsilon  method  yields  the  necessary 
conditions  of  optimality  obtained  by  applying  Pontryagin's 
minimum  principle  [Refs.  15  and  16]: 


x*(t)    =  |£  [x*(t),u*(t),t],  (2.5) 

~  dp        ~  ~ 


P*(t)    -  -  |£  [x*(t),u*(t),t],  (2.6) 

~  a  X        ~  ~ 


and 


H[x*(t),u*(t),p*(t),t]   <   H[x*(t),u(t),p*(t),t]  (2.7) 

•w  >W  "V  •**  •>•  **" 


where  H  is  the  Hamiltonian  function  defined  as 
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H[x(t),u(t),p(t),t]  =  g[x(t),u(t),t]  +  pT(t)f[x(t),u(t),t], 

(2.8) 

p(t)  is  the  costate  or  adjoint  vector,  u*(t)  is  an  extremal 
control  vector,  and  x*(t)  is  an  extremal  trajectory.   The 
assumptions  made  are  that  the  minimization  problem  has  a 
unique  solution  with  x(t)  absolutely  continuous  for  each  e, 
and  that  f  and  g  are  continuously  differentiable  in  x  and 
u  [Ref.  11].   Thus,  under  appropriate  assumptions,  it  can 
be  shown  that  as  e  ->  0,  the  epsilon  method  yields  the 
results  of  Pontryagin's  minimum  principle.   That  is,  if  the 
optimal  control  u*(t,e)  of  equation  (2.3)  exists  for  each 
e,  that  solution  will  approach  the  optimal  control  u*(t) 
of  equation  (2.2)  as  e  -*■   0. 

4.   State  Equation  Integration 

It  should  be  noted  that  the  epsilon  method  is  a 
non-dynamic  method  in  that  the  state  equations  are  not 
integrated  during  the  minimization  process.   Once  the  aug- 
mented performance  measure  has  been  minimized  a  check  on 
the  degree  of  satisfaction  of  the  state  equations  can  be 
obtained  by  integrating  the  state  equations  with  the  optimal 
control. 

v 

B.   MINIMIZING  THE  AUGMENTED  PERFORMANCE  MEASURE 
1.   Sequence  of  Unconstrained  Problems 

Once  the  augmented  performance  measure  (2.3)  is 
formulated,  any  unconstrained  optimization  algorithm  can 
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be  applied  to  It.   A  sequence  of  unconstrained  problems 

referred  to  as  sub-problems  is  solved.   In  each  sub-problem 

e  is  held  constant  and  a  minimization  is  performed  until 

some  stopping  criterion  is  satisfied.   At  this  point  e  is 

decreased  and  a  new  sub-problem  is  commenced  using  the 

optimum  trajectory  found  in  the  previous  sub-problem  as  a 

first  guess.   In  this  manner  a  sequence  of  sub-problems  is 

solved  until,  if  convergence  occurs,  some  overall  stopping 

criterion  is  satisfied. 

2.   Unknowns  and  Time  Discretization 

The  states  and  controls  can  be  approximated  by  any 

orthogonal  expansions.   The  coefficients  in  these  expansions, 

along  with  all  free  end  conditions,  become  the  parameters 

or  unknowns  in  the  optimization.   A  functional  expansion  of 

the  form  (2.9)  is  convenient  because  it  is  continuous  and 

the  period  can  be  selected  so  that  the  value  of  the  expansion 

is  zero  at  the  end  points.   Since  the  problems  solved  involve 

time-invariant  systems,  t  is  selected  as  zero  and  the  states 

o 

and  controls  are  written  as 


y(t)  = 


x(t) 
u(t) 


=  y 


Co  +  y(T)  -  y(0)   t  +  d 

T 


sin 


TTt 


T 
sin   2Trt 


T 


sin  Mut 


(2.9) 


where  M  is  the  number  of  harmonics  used  and 
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D 


D   = 


•  x 


D 
~u 


(2.10) 


is  an  (n  +  a)    x  M  matrix  of  coefficients.   The  derivative 
of  the  expansion  of  the  state  vector  given  in  equation  (2.9) 
with  respect  to  time  is  required  and  is  given  by 


x(t)  =  =- 


x(T)  -  x(0) 


T 


+  D 

~x 


7p  cos  -^p 


2tt  „„0  2ut 

ycos  -tjt- 


(2.11) 


The  objective  is  to  find  the  D  matrix  along  with  the  values 
of  the  free  end  conditions  which  minimize  equation  (2.3) 
for  a  given  e.   In  order  to  perform  this  minimization,  the 
time  interval  T  is  divided  into  (K  -  1)  sub-intervals  each 
of  duration  At  so  that  there  are  K  discrete  time  points. 
The  augmented  performance  measure  given  by  Equation  (2.3) 
is  written  as 

r  T  2       T  2 

J  (x,u,e)  -  - \  f   [x,(t)-f,(x,u)]  dt  +  f   [x5(t)-fP(x,u)]  dt 
a  ~  ~      ^L-'o        j.~~       j  q  a  <l   ~  ~ 


+  ..  .  + 


^[in(t)-fn(x5u)]2dt] 


(2.12). 


/'  gCx,u] 


dt  +  h[x(T),T]. 
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Suppressing  the  arguments  for  clarity,  equation  (2.12)  is 
expanded  to  yield 


.At 


KAt 


I/O    L     1  JtA  ^     1  ■'(K-DAt   X  X 


/-At  r2At  9  r     KAt 

/        (xp-f9)Mt  +/  (x„-fp)^dt  +  •••  +  /  (iL-fp) 

'0  d     d  JAt  d     2  -TK-lUt       *     2 


_dt 


(K-l)At 


/-At  _  /-2At  _ 

+     /        (x -f  )2dt  +/         (x . -f  ) 2dt  + 
/0  n     n  ^At  n     n 


(2.13) 


r    KAt 

y    ( 


Vfn)2dt] 


/-At  >-2At  /- 

/        gAt  +/  gAt  +   •••   +  / 

•'O  yAt  -TK 


KAt 


(K-l)At 


gAt  +  h[x(T),T] 


which  can  be  approximated  by 


a    ft 

+ 


f 


Xg 


1^(0)  -  ^WO),^)]!    ^+  |Xl(At)  -  yxCAtJ^CAt)]]    ^ 

2 

•  +/x1([K-l ]At)   -  f1[x([K-l]At),u((K-l)At)]j     ^  + 

0)  -  f2[x(0),u(0)]j    ^+  jxgCAt)  -  f2[x(At),u(At)]j     ^ 

•-+{xp((K-l)At)  -  fp[x((K-l)At),u((K-l)At)])    ^  +   •  ■  •   + 

l^  22~  "  2  (2'l4) 

0)  -  fn[x(0),u(0)]}     ^+  |xn(At)  -  fn[x(At),u(At)]}     ^ 

•  +  /^((K-DAt)  -  fn[x((K-l)At),u((K-l)At)]}     ^  + 
|g[x(0),u(0)]}At  +/g[x(At),u(At)]JAt  +   ...+ 
jg[x((K-l)At),u((K-l)At)]JAt  +  h[x(T),T]     .' 


fc 
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Hence,  the  augmented  performance  measure  can  be  written  as 


Ja  =  W;L2  +  w22  +  ...  +  wQ2  (2.15) 


wTw  (2.16) 


where  w  is  a  Q  x  1  column  vector.   The  first  K  elements  of 
w  are 

wk  =  [X;L[(k-l)At]  -  f1[x((k-l)At)1u((k-l)At)]}[^]%  , 

(2.17) 
k  =  1,2,..., X, 

etc.   The  form  of  equation  (2.16)  is  convenient  for  computer 
programming  the  epsilon  method  and  for  the  derivation  of 
the  minimization  techniques  that  follow.   In  minimum  time 
problems  where  the  performance  measure  is  given  by 


=  /  dt. 


(2.18) 


one  element  of  w  of  the  form 


wQ  =  C(K-l)At]%  (2.19) 


is  used  to  represent  equation  (2.18).   In  these  type  problems 
the  number  of  time  points  K  is  held  constant  during  the 
minimization  in  order  to  keep  the  dimensions  of  all  vectors 
and  matrices  constant  and  the  time  interval  At  is  minimized. 
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The  values  of  the  states  and  controls  required  in 
Equation  (2.1*J)  are  obtained  by  evaluating  equations  (2.9) 
and  (2.11)  at  each  time  point  t  =  (k  -  l)At  where  k  =  1,2, 
...,K.   Written  in  discrete  form  equation  (2.9)  is 


y[(K-l)At]-y(0) 

y[(k-l)At]  =  y(0)  +  - (k-1)  +  D 

K-l 


sin 

sin 


sin 


Tr(k-l) 
K-l 

2ir(k-l) 
K-l 


MTT(k-l) 

K-l 
(2.20) 


and  equation  (2.11)  is 


x[(k-l)At]  = 


x[(K-l)At]-x(0) 
(K-l)At 


+  D 
~x 


IT  TT(k-l) 

(K-l)At  COS 


2tt 
(K-l)At  C0S 


K-l 

2TT(k-l) 


K-l 


Mtt        MTt(k-l) 
(K-l)At  COS    K-l 


(2.21) 


A  vector  of  unknowns   c  is  formed  and  is  given  by 


c"  -  (d1}1>  dlj2>  •->    dljM>  d2,l'  d2,25  ■•"'  d2,M'  "' 


dn+Jt,l'  dn+A,2J  '••,  dn+Jl,M'  zl5  Z2>      ">    V  At) 

(2.22) 
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where  d±  .  is  the  element  in  the  i   row  and  jth  column 
of  D  and  z_,  z2,  ...,  zp  represent  P  free  end  conditions 
some  of  which  occur  at  t  =  0,  and  others  at  t  =  T.   Some 
of  the  z  ' s  correspond  to  states  and  others  to  controls. 
The  last  element,  At,  is  present  only  if  time  is  to  be 
minimized.   The  c  vector  consists  of  L  elements  where 


L  =  (n+Jt)  x  M  +  P  (2.23) 

for  all  problems  except  minimum  time  problems  and 

L  =  (n+£)  x  M  +  P  +  1  (2.24) 

for  minimum-time  problems. 

With  the  states  and  controls  given  by  equation 
(2.20)  and  the  augmented  performance  measure  given  by 
equation  (2.16),  the  problem  has  been  transformed  into  a 
parameter  optimization  problem  with  the  unknowns  given  by 
c  (2.22). 

3.   Minimization  Techniques 

The  methods  which  have  received  attention  in  the 
literature  for  finding  c*  which  minimizes  the  augmented 
performance  measure  given  in  equation  (2.3)  are  the  gradient 
method  and  a  "modified"  Newton-Raphson  method  (MNR) . 

The  gradient  method  has  been  investigated  by  J. 
Taylor  [Refs.  11  and  12]  and  L.  Taylor  [Ref.  8]  with 
unsatisfactory  results.   These  investigators  report  that 
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in  non-linear  problems  the  gradient  method  frequently 
obtains  false  minima  and  requires  considerable  computation 
time  compared  to  other  methods. 

An  MNR  method  in  which  certain  second-order  terms 
present  in  the  full  Newton-Raphson  (FNR)  method  are 
neglected  has  enjoyed  greater  success  and  requires  less 
computation  time  than  the  gradient  method  [Refs.  8,  11,  and 
12].   However,  in  Ref.  8  difficulties  are  reported  with  the 
MNR  method  in  non-linear  problems.   J  often  begins  an 
oscillation  after  two  or  three  iterations  and  does  not 
settle  to  a  minimum.   In  the  problems  solved  herein  the 
same  oscillations  have  been  observed  when  the  MNR  algorithm 
has  been  used.   Convergence  to  the  minimum,  when  it  does 
occur,  is  typically  very  slow.   Typical  performance  of  the 
MNR  method  is  shown  in  Figure  1. 
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Figure  1 


Augmented  performance  measure  vs.  iteration 
number-MNR  method 
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The  FNR  method  has  not  been  used  by  other  investi- 
gators with  the  epsilon  method  because  of  the  increased 
computer  time  for  each  iteration,  the  additional  storage 
space  required,  and  a  significantly  increased  analytic 
workload  involved  in  deriving  second  partial  derivatives. 
Because  of  the  poor  performance  of  the  MNR  method  on  problems 
of  the  type  solved  herein,  the  FNR  method  has  been  investi- 
gated in  detail  in  this  work.   With  careful  programming  the 
computer  time  for  each  iteration  and  storage  space  required 
for  the  FNR  method  has  been  reduced  to  an  extent  which  makes 
the  method  computationally  feasible. 

M.   The  Full  and  Modified  Newton-Raphson  Equations 

The  FNR  equations  for  finding  c*  are  derived  here 
in  a  manner  which  permits  the  MNR  equations  derived  in  the 
literature  [Refs.  8  and  11]  to  be  obtained  by  neglecting 
a  term. 

The  augmented  performance  measure  is  expanded  in  a 
Taylor  series  including  up  to  second-order  terms  and  is 
written  as 


J  (c+Ac)  2  J  (c)  .+  (VJ)Tic  +  ^(Ac)T(Vn2J  )(Ac)    (2.25) 


where 
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and 


V  J 
c  a 


3J. 


3cT 


(2.26) 


>\ 

c 

3C-3C2 

V2J       = 

c     a 

80236^^ 

3c22 

32J. 


3cL3c, 


32J. 


3cT  3Cp 


•     • 


•     • 


32J 


9CiaCL 


•      • 


32J. 


a 


3Cp3cT 


32J. 


3  c, 


(2.27) 


Solving  for  the  increment  of  J  ,  we  have 

3, 


AJ   =  J  (c+Ac)  -  J  (c) 


(2.28) 


s  (V  J  )/Ac  +  %(Ac)T(V„2J_)_(Ac).       (2.29) 


21 


Applying  the  necessary  condition  for  a  minimum,  we  have 


3(AJ  ) 

JTrAcf-  ■   (Va'c  +  (Vc  Ja}c  A2  =  2       <2'30) 


which  when  solved  for  Ac  yields 


Ac  =  -(V  2J  )  -1  (V  J  )  .  (2.3D 

c  a  c     c  a  c 


If  the  augmented  performance  measure  is  written  as 


A    T 
J   =  ww  (2.14) 

3.  -*  *•* 


where  w  is  a  Q-dimensional  column  vector,  then 


Va  =  V^ 


=   2(Vcw)Tw 


and 


Vc2ja  =   2Vc[(Vc^)T^] 


=  2[(Vcw)T(Vcw)  +  (Vc2w)Twl 


The  matrix  V  w  is  given  by 
c~ 


(2.32) 


(2.33) 


22 


3w. 


aw. 


3w, 


3c 


3c, 


3cn 


Vc? 


3w, 
3cl 


3w, 
3cT 


3w, 


(2.34) 


3w 
jTcT. 


Q 


3w 
3~cT 


Q 


3w 

3c; 


Q 


V     w   is   a  three-dimensional   array   composed  of  L  matrices   each 

C         ~  Q 

32wi 


of  dimension  Q  x  L  which  has  as  its  ijk  th  element 


that  is 


c  ~ 


"w  = 


3  w- 


3c. 


32w, 


3c 


3  w 


Q 


3c 


32w, 


2 

3  w. 


3c13c2 


3c13cL 


32w, 


32w, 


3c13c2 


3c,3cL 


32w 


Q 


32w 


Q 


3c-,3c2 


3c-,3cL 


32w 


Q 


3  w 


Q 


SCpSc, 


3c. 


2 

3  w 


Q 


3cr 3c. 


32w. 


3c23cL 
32w. 


3c23cL 


32w 


Q 


3c23cL 


32w 


SL 


3cL3c2 


3c,  3c.  ' 


3c, 


(2.35) 
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Substituting  equations  (2.32)  and  (2.33)  into  equation  (2. 3D, 
we  have 

As  ^Y^Y's +  cV~)s\^1(ve?)2T  »c  •   (2-36) 

This  is  the  full  Newton-Raphson  equation.   The  modified 
Newton-Raphson  equation  can  be  obtained  by  neglecting  the 
second  term  in  the  inverse  in  equation  (2.36),  which  yields 

As  -  -  f(Y)cT(vc^)£j"1(Y)sT  ~s  (2,37) 

Several  comments  concerning  equations  (2.36)  and  (2.37)  are 
in  order: 

a.  the  term  V  w  given  by  equation  (2.3*0  is  a  Q  x  L  matrix; 

b.  the  term  (V  w)  (V  w)  is  a  symmetric  L  x  L  matrix; 

m 

c.  the  term  (v  w)  w  is  an  L  x  1  vector; 

d.  the  term  v  2w  given  by  equation  (2.35)  is  a  Q  x  L  x  L 

c  ^ 

three-dimensional  array; 

2  T 

e.  the  term  (V  w)  w  is  a  symmetric  L  x  L  matrix; 

?   T     "" 

f.  (V  w)   is  defined  as  the  three-dimensional  array 

obtained  by  transposing  each  individual  Q  x  L  matrix 

given  in  equation  (2.35); 

2   T 

g.  the  result  of  the  operation  (V  w)  w  is  defined  as 

an  L  x  L  matrix  in  which  the  i   column  is  the 
product   I  r—  (V  w)    w  . 


n —  (v  w)T 

[   3ci    c~J J 
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5.   A  Scalar  Illustration  of  the  MNR  and  FNR  Methods 
The  potential  importance  of  the  second-order  term 
neglected  in  the  MNR  equations  can  be  illustrated  with  a 
simple  scalar  problem  in  function  minimization.   The  Newton- 
Raphson  equation  to  minimize  a  function  f(x)  takes  the  well 
known  form 


Ax  =  -  f^fj-   .  (2.38) 


If 


f(x)  =  w2(x)  (2.393 


which  is  the  form  of  equation  (2.16)  with  w  taken  as  a 
scalar,  then 

f'(x)  =  2w(x)  g  (x)  (2.40) 

and 

.  f"(x)  =  2  [g  (x)]2  +  2w(x)  S-5  (x)   .     (2.41) 


dx 


The  FNR  equation  (2.36)  is 


-w(x)  g  (x) 

Ax  = g (2,1<2) 

d  w 


[  S  <*>]  +  »<*>  t?  (x) 


whereas  the  MNR  equation  (2.37)  is 
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-»<*>  &  <*> 


Ax  =  (2.13) 


"(x)  (2  214) 

dx  vx; 


Applying  equations  (2.42)  and  (2.44)  to  the  function 


f.(x)  =  (x-1)4  +  1  (2.45) 


in  which 


w(x)  =  [(x-1)4  +  1]%   ,  (2.46) 


we  obtain 

Av  =  _ 

3 


Ax  =  -  |  (x-1)  (2.47) 


for  the  FNR  algorithm  and 


Ax  =  -  &=^ — ±-i  (2.48) 

2(x-l)3 


for  the  MNR  algorithm.   Tables  1  and  2  show  the  first  few 
iterations  by  both  methods  from  an  initial  guess  of  x(0)  =  3 
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MNR  Equation  (2.48) 

Iteration 

X 

Ax 

f(x) 

1 
2 

3 

4 

3.000 

1.937 

0.862 

190.932 

-1.062 

-1.076 

190.070 

145.121 

17.000 
1.772 
1.000 
1.329  x 

109 

Table  1 

FNR  Equat: 

Lon  (2.47) 

Iteration 

X 

Ax 

f(x) 

1 

3.000 

-0.667 

17.000 

2 

2.333 

-0.444 

4.160 

3 

1.889 

-0.296 

1.625 

4 

1.593 

-0.197 

1.124 

5 

1.395 

-0.132 

1.024 

6 

1.263 

-0.088 

1.005 

7 

1.175 

-0.058 

1.001 

8 

1.117 

-0.039 

1.000 

9 

1.078 

-0.026 

1.000 

Table  2 


Clearly  the  MNR  equation  causes  x  to  diverge  after  an  initial 
period  of  convergence  while  the  FNR  equation  causes  x  to 
approach  the  minimum. 

6.   Computation  Experience  With  the  FNR  and  MNR  Methods 
The  performance  of  the  MNR  method  in  the  preceding 
scalar  example  is  typical  of  the  performance  observed  by  the 
author  in  large  problems.   However,  the  FNR  equation  is  also 
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not  uniformly  effective  when  used  exclusively  in  large 
problems.   Fortunately,  the  areas  of  effectiveness  of  the 
two  methods  are  complementary. 

In  order  to  discuss  the  effectiveness  of  the  two 
methods  it  is  convenient  to  define  two  areas  in  the  minimiza- 
tion process.   Initial  behavior  refers  to  the  behavior  of  J 

a 

during  the  first  two  or  three  iterations  in  a  sub-problem. 

Terminal  behavior  refers  to  the  behavior  of  J  after  the 

a 

first  two  or  three  iterations  within  the  same  sub-problem. 
The  following  behavior  has  been  observed. 

a.  Initial  behavior:   The  MNR  equation  outperforms 
the  FNR  equation  in  this  area.   The  ability  of  the  FNR 
equation  to  minimize  J  is  very  sensitive  to  the  starting 
value  of  the  unknowns  (c).   With  the  values  of  c  far 
removed  from  the  optimum,  the  FNR  equation  generally  causes 
J_  to  increase  rapidly  and  diverge  from  the  minimum.   The 

cl 

MNR  equation  on  the  other  hand  is  relatively  insensitive  to 

the  starting  c  and  can  usually  be  counted  on  to  move  J 

~  a. 

toward  the  minimum  for  at  least  one  or  two  iterations. 

b.  Terminal  behavior:   As  the  minimum  is  approached, 
the  MNR  equation  produces  the  behavior  shown  in  Figure  1. 
The  FNR  equation,  however,  generally  becomes  extremely 
effective  in  rapidly  finding  the  minimum. 

7.   A  Combination  FNR-MNR  Minimization  Method 

The  obvious  approach  suggested  by  the  previous  obser- 
vations is  to  devise  an  algorithm  which  minimizes  by  the  MNR 
equation  initially  in  a  given  sub-problem  and  switches  to 


28  i   N 


the  FNR  equation  at  some  appropriate  point  in  the  iteration 
process.   Such  an  algorithm  is  presented  in  Section  IV. 
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III.   AN  INEQUALITY  CONSTRAINT  PENALTY  FUNCTIONAL 

In  this  section  a  new  penalty  functional  is  introduced 
for  state  and  control  inequality  constraints  of  the  form 


x,   <  x  (t)  <  x,    ,   1=1, 2,..., I  <  n  ,   t  e  [to,T]  ,  (3.1) 


u,   <  u,(t)  <  u    ,   j=l,2,...,I  <  I   ,   t  e  [t  ,T]  .  (3.2) 


All  state  and  control  inequality  constraints  encountered  in 
the  problems  solved  herein  are  of  this  type.   The  difficul- 
ties encountered  with  existing  penalty  methods  which  led  to 
the  use  of  a  new  functional  are  related.   The  new  penalty 
functional  has  performed  well  in  computation  and  is  used 
exclusively  in  the  solution  of  the  problems  presented. 
Additionally,  several  desirable  theoretical  properties  of 
the  proposed  penalty  functional  are  presented. 

A.   INTERIOR  PENALTY  METHODS 
1.   Past  Research 

In  Ref.  10  Jones  and  McCormick  present  a  number  of 
theoretical  results  concerning  interior  penalty  functionals 
of  the  Fiacco-McCormick  type  [Ref.  9]  in  conjunction  with 
the  epsilon  method.   If,  for  example,  a  state  or  control, 
denoted  by  y(t)  for  generality,  is  constrained  by 

y(t)  <  Y   ,    t  e  [t  ,T]  ,  (3.3) 
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a  Fiacco-McCormick  penalty   functional   [Ref.    9]   of  the    form 


/ 


T 


1  -  y(t) 


dt 


(3.H) 


is  added  to  the  augmented  performance  measure.   The  behavior 
of  the  integrand  of  expression  (3.^0  for  a  fixed  time 
t  e [t  ,T]  as  the  positive  weighting  factor  r  approaches  0 
is  shown  in  Figure  2. 


Figure  2 

Fiacco-McCormick  penalty  function  vs.  constrained  variable 

for  a  fixed  time 


31 


If  penalty  functionals  of  this  type  are  added  to  the  perform- 
ance measure,  it  can  be  shown  [Ref.  10]  that  as  r  approaches 
0  and  e  approaches  0,  the  epsllon  method  yields  Pontryagin's 
minimum  principle.   The  development  parallels  and  augments 
Balakrishnan's  work  [Refs.  13  and  14]  without  inequality 
constraints.   No  computational  results,  however,  are 
presented. 

2 .   Computational  Experience 

A  simple  problem  involving  one  state  variable  and 
one  constrained  control  was  attempted  using  the  epsilon 
method  and  a  penalty  functional  of  the  form  given  by 
equation  (3-^)-   The  optimal  control  was  on  the  constraint 
boundary.   The  algorithm  was  unable  to  solve  the  problem 
by  either  the  FNR  or  MNR  method  from  a  variety  of  starting 
points.   Once  the  control  penetrated  the  constraint  boundary 
for  a  finite  time  interval,  the  algorithm  failed  on  the  next 
iteration.   The  value  of  r  required  to  keep  the  control 
admissible  for  all  t  e  [t  ,T]  throughout  the  iteration 
process  was  large,  resulting  in  the  augmented  performance 
measure  being  dominated  by  the  Fiacco-McCormick  penalty 
term.   As  a  result,  the  optimal  solution  [u*(t,e,r)]  to  the 
augmented  problem  could  not  be  made  to  approach  the  optimal 
solution  [u*(t)]. 

B.   A  NEW  PENALTY  FUNCTIONAL:   COMPUTATIONAL  PROPERTIES 
1.   The  Form  of  the  New  Penalty  Functional 

Consider  a  control  or  state  y(t)  which  is  subject  to 
a  constraint  of  the  form 
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y(t)  e  [yL,yM]   ,    t  e  [tQ,T]   .        (3.5) 


A  penalty  functional  of  the  form 


t  r2y(t)  -  yM  -  yTi2K 


-T  r2y  t)  -  yM 


dt    (3.6) 


where  K  Is  a  positive  integer  is  added  to  the  augmented 
performance  measure.   The  effect  of  K  can  be  seen  from 
Figure  3  which  shows  the  integrand  of  equation  (3-6)  as  a 
function  of  y(t)  for  a  fixed  time  t  e  [t  ,T]. 


increasing  K 


Figure  3 
New  penalty  function  vs.  constrained  variable  for  a  fixed  time 
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A  functional  of  the  form  given  by  equation  (3.6)  is  added 
to  the  augmented  performance  measure  for  each  inequality 
constraint  of  the  form  given  by  equations  (3.1)  and  (3'.2). 
For  I   control  constraints  and  I  state  constraints  of  this 

B  S 

form  the  total  augmented  performance  measure  for  the  epsilon 

method  written  for  time  invariant  problems  with  t  =  0  is 

o 

T 
Ja  =  J+/     |  ||i-  f(x,u)||2dt 

(3.7) 


-TV  ^c  ,-On    /'+->_,-,,   _,,     -,r>v  ' 


Jo  l. 1-iL   ui  -uiT  J      M   xi.-xi    J   / 


JM         JL 


1=lL    \'\ 


dt 


J  +  I  Js   +  rJp  (3.8) 


The  ".two-sided"  feature  of  the  penalty  functional  makes  it 

especially  suited  to  constraints  of  the  form  given  by 

equations  (3.1)  and  (3-2).   In  effect,  two  inequality 

constraints  are  included  in  one  penalty  term. 

2.   Computational  Strategy 

The  power  K  is  increased  gradually  in  numerical 

computation  in  the  same  manner  as  e  is  decreased.   Thus, 

increasingly  refined  boundaries  to  the  admissible  region 

are  provided.   Both  e  and  K  are  held  constant  within  a 

P 

sub-problem  and  are  altered  between  sub-problems.   The 

weighting  factor  r,  which  is  required  to  provide  an  overall 

weighting  among  J,  J  ,  and  J  ,  is  held  constant  throughout 

s       p 

the  entire  problem. 

3^ 


Computational  results  with  this  penalty  functional 
have  been  excellent.   Up  to  eight  inequality  constraints  have 
been  treated  successfully  in  one  problem. 

C.   THEORETICAL  PROPERTIES  OP  THE  NEW  PENALTY  FUNCTIONAL 
1.   Introduction 

Three  desirable  properties  of  the  new  penalty 
functional  are  presented  here.   These  properties  and  their 
importance  are  discussed  followed  by  a  proof  of  each 
property. 

a.   Penalty  functionals  of  the  form  of  equation  (3.6) 
are  convex  on  Rn  where  c   e  Rn  is  defined  by 

cyT  =  Ca1,a2,...,aM,y(t0),y(T)]  (3-9) 


and 


y(T)-y(t  )        M      imr(t-t  ) 
y(t)  =  y(t  )  +  =  .  °  (t-t  )  +  Vam  sin   _  .  °       (3-10) 


T-t     v   o'   Z_-»  m      T-t 
o  m=j 


It  is  desirable  that  the  augmented  performance  measure  be 
convex  in  the  unknowns  of  the  minimization  to  insure  that  a 
global  minimum  is  attained.   If  it  can  be  shown  that  the 
inequality  constraint  penalty  functionals  (3-6)  are  convex, 
then  the  addition  of  any  number  of  these  functionals  does 
not  destroy  a  convexity  condition  which  exists  without  these 
terms,  because  the  sum  of  convex  functionals  is  also  convex. 
Indeed,  the  addition  of  terms  of  the  type  given  in  equation 
(3.6)  may  create  a  convexity  condition  where  one  does  not 
exist  without  the  terms. 
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b.   If  for  a  fixed  c  given  by  equation  (3-9),  the 
expansion  (3.10)  of  a  constrained  state  or  control  is  inad- 
missible by  a  finite  amount  e  (e   >  0)  at  txe(t  ,T),  its 
associated  penalty  functional  (3.6)  is  unbounded  as  K  •*■«*. 
This  means  that  as  K  ■+  »,  the  contribution  of  a  penalty 
term  (3 .6)  to  the  augmented  performance  measure  for  an 
inadmissible  state  or  control  becomes  very  large.   Therefore 


if  J„  is  being  minimized  under  the  condition  of  ever  increas- 

el 

ing  K  ,  the  constrained  states  or  controls  must  at  least 
approach  admissibility. 

c.   If  for  a  fixed  c  given  by  equation  (3-9),  the 
expansion  (3.10)  of  a  constrained  state  or  control  lies 
completely  within  the  admissible  region,  its  associated 
penalty  functional  (3-6)  has  limit  zero  as  K  ■*■  °°.   The 
significance  of  this  result  is  that  penalty  terms  of  the 
form  given  by  equation  (3-6)  will  add  less  and  less  to  the 
augmented  performance  measure  as  K  is  increased  for  states 
and  controls  that  are  completely  admissible. 
2.   Convexity 

Property  a  discussed  above  is  shown  here.   The 
theorem  to  be  proved  follows. 

Theorem  1.   If  a  constrained  state  or  control  y(t) 
is  bounded  for  t  e  [t  ,T]  and  is  given  by 


y(T)-y(t  )         M      rmr(t-t  ) 
y'1  !  -  "   !    nx_+  ~  Ct-tJ  +  £am  sto   T_t        (3.10) 


;  )  + a  z  °  (t-t  )  +    Ta    sin  — 7=-r- 

0      T-tQ       o    ^      T-tc 


where  c   e  Rn  is  defined  by 

~y 
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cyT  =  [a1,  a2,  ...,  aM,  y(tQ),  y(T)]  ft   0  ,    (3-9) 


then  the  penalty  functional 


Tr2y(t)-yM-yLi 


Jjy,Kj    =   r  f     — 

o 


p     p       t{  L    yM-yL 


2K 

p  dt  (3.6) 


is  convex  on  R  .   The  constants  y„  and  y  define  the 
admissible  region  for  y(t)  and  r  is  a  constant. 


Proof.   Consider  the  case  of  K  =1.   Equation  (3.6) 
becomes 


J  (y,K  )  =  £ p  I    [2y(t)-(y  +y  )]2  dt.  (3-11) 

p    p    (v  -v  )  tJ 

vyM  yI/    o 


Let 


a  yM  +  yT 

d  k     M  -  L  (3.12) 


and 


A        ^r  ,,   HON 

r  =  o  •  (3-13) 

(yM-yL} 


Substituting  equations  (3-12)  and  (3-13)  into  equation  (3.11), 
we  obtain 
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1 

=   ro     /   Cy(t)    -   d]I 


dt 


(3.14) 


T 

-  rn      f  [y2(t)   -  2y(t)d  +  d2]   dt    .  (3.15) 

Substituting  the   expansion   (3.10)    into   equation    (3.15),   we 
obtain 


J    =  r 
P        o 


/-V  y(T)-y(t  )  M  imr(t-t  )-,2 

/  {ro)  +        TV     <**o>  +  E   a^  sin  -^ 


m=l 


(3.16) 


r  y(T)-y(t  )  M  iror(t-t  ) 

* kV  +        T-t  «Hb0)  +   E   a^  sin  -^-5- 

I  o  m=l  o 


d  +  d2l 


dt. 


Equation  (3.16)  may  be  written  as  a  quadratic  functional 
of  the  form 

T 

Jp    =    r0       /    [i<V~l(t)Sy>    +   <2y>?(t»    +    B]  dt     (3'17-) 


n 


where   a(t)  e  R     is 


a(t)T  = 


-2d 


sin 


ir(t-tQ)       sin2Tr(t-t0) 


T-t       *         T-t 
o  o 


sinMTT(t-t  ) 


T-t 


t-t      t-t 

(1 °)  ° 

,u  T_t  ;,T_t 

o  o  -> 


(3.18) 
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Q1(t)  is  the  outer  product  given  by 


Q_(t)  =  silMl) 
-1  Hd2 


T 


(3.19) 


and 


g  =  d< 


(3.20) 


At  an  arbitrary  fixed  time  t*e[t  ,T],  the  integrand  of 
equation  (3.14)  is 


i2      1 


[y(t*)-d]  =  2  <2y»Si(t*)£y>  +  <Sy^(fc*)>  +  P  •    (3.21) 


The  first  term  on  the  right  side  of  equation  (3.21)  is 


iKvV^Sy) =  [y(t< 


y(T)-y(t)         H      nnr(t#-t)n2 
J  +  -ktt— 2-  (t»-t0)  +  £  V^-^T-2- 


T-t 


m=l 


(3-22) 


Since  the  terms  in  the  finite  expansion  of  y(t#)  given  in 
equation  (3.10)  are  linearly  independent  and  analytic,  y(t#) 
is  different  from  zero  almost  everywhere  and 


y  (t#)  >  0  ,    t*  e  [tQ,T],    cy  ?   0, 


(3.23) 


Thus, 


h  <£  .Q^t^c  >  >  0,  t«  E  [tQ,T],   cy  ft   0, 


(3-2*0 
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and  Q-j^Ctjj)  is,  therefore,  positive  semi-definite,  at  least, 

and  positive  definite  almost  everywhere  for  any  c  ^0. 

~y 

Applying  Theorem  4.5  of  Reference  17  (p.  2?)1,  the  function 

given  in  equation  (3-21)  is  convex  on  Rn  at  t  =  tr 

Next,  consider  the  case  where  K  is  any  positive 

P 

integer.   In  Appendix  E  the  following  theorem  is  proved: 
if  f(x)  is  convex  on  Rn  where  x  e  Rn  and  f(x)  >_  0,  then 
r  (x)  is  convex  on  R  where  K  is  any  positive  integer. 
Since 


[y(t»)  -  d]2  >  0,  (3.25) 


is  convex,  it  follows  immediately  that 


2K 
[y(t,)  -  d]  P  (3.26) 


is  convex  on  R  at  a  fixed  time  tx  e  [t  ,T].   Since  y(t)  is 
bounded  by  assumption  for  t  e  [t  ,T],  it  follows  that  for 


Let  f  be  a  twice  continuously  differentiable  real- 
valued  function  on  an  open  convex  set  c  in  R  .   Then  f  is 
convex  on  c  if  and  only  if  its  Hessian  matrix 

2 
Qx  =  (q^Cx)).  q1;,<x)  =  f^-  tev    ...,  Kn) 

is   positive    semi-definite   for  every   x   e   c.      A  quadratic   function 

f(x)    =   %  <x,Qx>     +     <x,a>     +   a 

where  Q  is  a  symmetric  n  x  n  matrix,  is  convex  on  R  if  and 
only  if  Q  is  positive  semi-definite. 

l\0 


any  finite  positive  integer  K  ,  the  expression  (3.26)  is 
bounded  for  tx  e  [t  ,T].   By  Theorem  4,  (p.  536)  of 
Reference  18 


T  2K 

f   [y(t)  -  d]  p  dt  (3.27) 


is  convex  on  R  .   Hence  equation  (3.6)  is  convex  on  Rn  for 

all  finite  positive  integer  values  of  K  . 

P 

3.   Behavior  of  the  New  Penalty  Functional  for  an 
Inadmissible  Constrained  State  or  Control 

Property  b  is  shown  below. 


Theorem  2.   Assume  y(t)  is  bounded  and  given  by 
equation  (3.10)  for  t  e  [t  ,T]  where  c  ¥   0,  as  defined  by 
equation  (3-9).   If  for  a  given  y(t) 


y(t*)  >  yM  +  ep  (3-28) 


or 


y(t*)  <  yL  -  eP  (3,29) 


■'"Let 


If(x)  =   /f(t,x(t))  dt. 


Let  T  be  of  finite  measure.   Let  f(t,x)  be  a  finite  convex 
function  of  x  for  each  t  and  a  bounded  measurable  function 
of  t  for  eachMx.   Then  I_  is  a  well-defined  finite  convex 
function  on  L°°(T)  which  Is  everywhere  continuous  with 
respect  to  the1  uniform  norm. 


HI 


at  some  time  t#  e  (t  ,T),  where  e   >  0,  then 


Limit  J  (y,K  )  =  ~  (3-30) 

K  -**>        p    p 
P 


where 

T  r2y(t)-y  -y  -,2K 
J  (y,K  )  =  r  f       ?$—k 

o 


p 


dt.       (3-6) 


In  order  to  prove  this  theorem  the  following  lemma  is 
required. 

Lemma  1 .   Assume  y(t)  is  bounded  by 


-<*>  <  M1  <  y(t)  <  M2  <  oo  (3-31) 


and  is  given  by  equation  (3.10)  for  t  e  [t  ,T]  where 

c  ^  0  as  defined  by  equation  (3-9).    Then,  if 
~y   ~ 


y(t»)  >  yM  +  ep  (3-28) 


at  some  time  t«  e  (t  ,T)  for  any  e_  >  0,  there  exists  a 
*     o  p 

6  >  0  such  that 


y(t)  >  yM  +  Z-f  (3.32) 


for  the  finite  time  interval 


t,  -  «  <  t  <  t,  H.  (3-33) 


H2 


•■ 


Similarly,  if 


y(t»)  <  yL  -  £p 


(3.29) 


at  some  time  tx  e  (t  ,T)  for  any  e   >  0,  there  exists  a 
6  >  0  such  that 


y(t)  <  yL  -  -£■ 


(3.31) 


for  the  finite  time  interval 


tg  -  6  £  t  £  t#  "+  6 


(3.33) 


Proof  of  the  lemma.  First,  it  is  necessary  to  show 
that  the  coefficients  in  equation  (3.10)  are  bounded;  that 
is 


|aj  <  M3,   m  =  1,2,.  . .  ,M 


(3.35) 


where   NU   >   0 .      To  this   end  consider 


M 


/?  J?r  y(T)-y(t  ) 


a  sin 


imKt-t  )-i2 
o 


m        T-t 
m=l  o 


dt 


(3-36) 


=    <SyS22y> 


(3-37) 


^3 


l 


where  c   is  given  by  the  definition  (3-9).   Q?  is  a 
symmetric  matrix  given  by 


Q2 


T-t 


2(T-to) 


IT 


2(T-tQ) 


0    . 


T-t 


o    .       .      0 


•  • 


2(T-tQ) 

2(T-tQ) 
2tt 


T-t 


2(T-tQ) 

Mi 

Mtt 


2 (T-t    ) 
o 

u 

2(T-tQ) 
2rf 


2(T-tQ) 

MlT 

Mo 

3 

*"*<> 

2 (T-t    ) 
o 

ir 

.  2(T-tQ) 

2tt 

• 
• 
• 

h2(T-tQ) 

Mtt 

Mo 

3 

Mo 

;3 

(3.38) 

where  the  terms  with  1   are  positive  if  M  is  odd  and  negative 
if  M  is  even.   Since  y(t)  is  bounded  by  assumption  and  the 
expansion  (3.10)  is  the  sum  of  M  +  2  linearly  independent 
terms,  we  have 


T 
0  <  f   [y(t)]2dt  <   M4 
t_ 


(3.39) 


where  Mm    >   0.      Using  equation    (3-37)    and  inequality    (3-39), 
we   obtain 


0    <     <2y'S2Sy>     -  M4 


(3.40) 
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Q2  is,  therefore,  positive  definite.   Using  Theorem  2.5 
of  Reference  15  (p.  52)  ,  we  have 


<  VS2£y>  -"  "Sy"2  C3.4D 


where  X_  >  0  is  the  smallest  eigenvalue  of  Q?  and 
ll£y||  =  ^/<c  ,c  >   .   Therefore,  using  the  inequality 
(3.^0),  we  obtain 


HSyll   1  T  *  (3'1,2) 


Since  a  ,   m  =  1,2,...,M  is  a  subset  of  c  ,  it  follows  that 


laml  1  M3>   m=l,2,...,M  (3.35) 

where  M  >  0. 

Now  consider  the  derivative  of  equation  (3.10) 
which  is 


y(T)-y(t  )    M  miT(t-t  ) 

*<*>  =    T-t     +  2>m  T?f  cos  -T3F-2-  •      (3.43) 
o      m=l      o         o 


Let  Q  =  (q^j)  be  a  symmetric  n  x  n  matrix.   Then  Q  is 
positive~definlte  if  and  only  if  there  is  a  k  >  0  such  that 

<v,Qy>  >  k||v||2 

for  all  v  in  Rn,  where  ||v||  =   «/<v,v>   is  the  Euclidean 
norm  of  v. 
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Taking  the  absolute  value  of  both  sides  of  equation  (3.43) 
and  applying  inequality  laws,  we  obtain 

y(T)-y(t  )      M     m^  imr(t-t  ) 

|y(t)  |  <  | T_t  °  1  +  |  E  am  g£-   cos      S  1    (3.44) 

o        m=l      o         o 

which  further  simplifies  to 

y(T)-y(t  )     M 

|yCt)|  1  1   T.t  °  I  +  E  Iaml  5^-  •      (3.45) 

o       m=l        o 

Applying  the  inequality  (3. 35),  we  obtain 

y(T)-y(t  )    M,tt    M 

|y(t)|  <  | T_t    1  +  tV   Em  (3.46) 

o  o  m=l 


<  MR  (3.47) 


for  all  t  e  (t  ,T)  where  IVU  >  0.  The  first  part  of  the 
lemma  as  expressed  by  the  inequality  (3-32)  will  now  be 
shown.   Let 


6  ^  -g—  (3.48) 

5 


as  shown  in  Figure  4.   Consider 

y(t)  =  y(t*)  +  J    y(x)  dx   .  (3-49) 

t » 
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y 

y(t) 

\  4 ! 

1 

/             / 

P 

/ 

/slope  =  M5           / 

6 

—      ^ 

\          \~"slope 

t 

* 

t 

Figure  4 
Constrained  Variable  vs.  time 


Applying  inequality  (3.^7),  we  have 


y(t)  >  y(t#)  -  M   |t,  -  t|  . 


(3-50) 


Consider  t  in  the  interval  t#  -  <S  <   t  <_  tx  +  6.   In  this 
interval  6  satisfies 


6  >  |t,  -  t|  . 


(3-51) 


Applying  inequality  (3-50),  we  have 


y(t)  >  y(t#)  -  M5  6 


(3-52) 
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Using  the  definition  (3.48),  we  obtain 


y(t)  >  y(t#)  -  -g.   . 


(3-53) 


Applying  inequality  (3.28),  we  obtain 


y(t)  >  yM  +  ep  -  -f 


(3.54) 


or 


yet)   >  yM  +  -#  , 


M 


(3.32) 


thus  proving  the  first  portion  of  Lemma  1.   The  second  portion 
of  the  lemma  as  given  by  inequality  (3-3*0  can  be  proved  in 
a  similar  manner.   It  is  possible  at  this  point  to  return 
to  the  proof  of  Theorem  2. 

Proof  of  Theorem  2.   If  inequality  (3-28)  applies, 
then  by  Lemma  1  inequality  (3*32)  is  true.   Considering  the 
integrand  of  equation  (3-6),  since  yM  >  y, ,  it  follows  that 


-2y(t)-yM-yLi2K 


yM  yL 


2(yM+  2  )  yM-yL~i 

yM"yL 


2K 


(3.55) 


for  t.  -  i  <  t  <  t|  +  i.   Further  simplification  yields 


H8 


<   \ 


2K   rer,+yM-y1 


F2y(t)-yM-yLn2Kp  >  r^ 
L   yM~yL   J     "  L  yM"y 


2K 


(3.56 


[ 


>  1  + 


yM  yLJ 


2K 


(3.57) 


Since  the  integrand  of  equation  (3.6)  is  nonnegative  for 
t  e  (t  »T)  and  r  >  0,  it  follows  that 


Vy-V  =  rft    [ 


t  r2y(t)-yM-yL 


O 


yM  yL 


2K 


dt 


(3.6) 


t/-d 


**'  r2y(t)-yM-yL 


yM  yL 


2K 


dt  . 


(3-58) 


Applying  inequality  (3-57),  we  have 


t»+6 


Vy'V   1  r  /  I"1  +  y~^ 

P    P      t/-6   L     yM 


LJ 


2K 


dt  . 


(3.59) 


The  integration  of  inequality  (3-59)  yields 


JpCy 


'V  2  r  [ 


e„  H2K 


1  + 


yM  yL- 


26 


(3.60) 


Since   y„   >   yT    and   e      >   0,    it   follows   that 
ML  p 


Limit 

K     -»•«> 

P 


r 

£~     1 

1 

+ 

P 

yM~yL-l 

n2K 


P  =      CO 


(3.61) 
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and  in  view  of  inequality  (3.60) 


Limit  _  /   „  v 

Kp+«  Jp(y>Kp)  =  °°  * 


(3-62) 


If  inequality  (3-29)  applies,  then  by  Lemma  1, 
inequality  (3-31*)  is  true.   Since  yM  >  y  ,  it  follows  that 


F 


2y(t)-yM-yTi    r2(yT-  -£)-y 


'M  'L 


yM  yL 


'L~  T;~JM"yL 
yM~yL 


(3.63) 


for  t*-<5<t<tx   +   6.      Simplifying,    we   obtain 


2y(t)-yM-yLl 

yM"yL 


<    - 


1   + 


yM  yLJ 


(3.64) 


Squaring  both  sides  of  inequality  (3.64),  we  have 


2y(t)-yM-yL 
yM~yL 


-.2 


-i2 


1  + 


yM-yL 


(3.65) 


Raising  inequality  (3-65)  to  the  K  th  power,  we  obtain 


2y(t)-yM-yL 
yM~yL 


2K 


>-[ 


e   n2K 


1  + 


yM  yLJ 


(3-57) 


which  is  identical  to  inequality  (3-57).   The  remainder  of 
the  proof  follows  inequality  (3. 57)  to  equation  (3.62) 
exactly.   The  theorem  is  proved  for  the  open  interval  t  e  (t  ,T) . 
The  extension  of  the  theorem  to  cover  the  closed  interval 
t  e  [t  ,T]  is  not  difficult. 
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4.   Behavior  of  the  New  Penalty  Functional  for  an 
Admissible  Constrained  State  or  Control 

Property  c  above  is  shown  below. 

Theorem  3-   Assume  y(t)  is  given  by  equation  (3.10) 

for  t  e  [t  ,T],   If  for  a  given  y(t) 

yL  +  eq  1  y(t)  1  yM  "  eq  (3'66) 


for  all  t  e  [t  ,T]  where  e   >  0,  then 


v-"  vy'V  *  °  •  (3-67) 


Proof.   From  inequality  (3-66)  it  can  be  seen  that 


yL  +  £q  1  yM  "  £q  *  (3'68) 


Since  y„  >  y,,  inequality  (3-68)  may  be  rarranged  to  the  form 


2en 
1  -  — g-  >   0   .  (3.69) 


The  inequality  (3-69)  will  become  useful  shortly. 

Starting  with  inequality  (3-66),  multiplying  through 
by  2,  and  subtracting  yM  and  yL,  we  obtain 


2<yL+Eq)-yM-yL  ±  2y(t)-yM-yL  ±  ^n-V"7*"*!,   •        (3-70) 
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Dividing  through  by  the  positive  quantity  yM  -  yL,  we  obtain 


2<yL+e  )-yM-yL   2y(t)-yM-yL    2(yM-e  )-yM-yL 


yM  yL 


yM  yL 


yM~yL 


(3.7D 


This  inequality  can  be  reduced  to 


1 9- 


yM  y 


M  JL 


r2y(t)-yM-yLn 

yM"yL 


2e. 


<  1  - 


yM  yL 


(3-72) 


In  view  of  inequality  (3.69),  inequality  (3.72)  may  be 
rewritten  as 


2y(t)-yM-yL 
yM"yL 


2e 


<  1  - 


yM  yL 


(3-73) 


Raising  inequality  (3-73)  to  the  2K  power,  we  have 


2y(t)-yM-yL 
yM"yL 


2K 


2e   H2K 


1  - 


yM"yL 


(3.74) 


Observing  that 


i[ 


2y(t)-yM-yL- 


yM  yL 


2K 


P  = 


2y(t)-yM-yL- 
yM"yL 


2K 


(3-75) 


we  have  from  inequality  (3.7*0 


2y(t)-yM-yL 
yM"yL 


2K 


2e 


1  - 


yM  yL-J 


2K 


(3.76) 


for  all  t  e  [t  ,T].   In  view  of  inequality  (3-76),  it  is 
easily  seen  that 
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Tr 

/[ 


2y(t)-yM-yL-,2K 
yM~yL 


/-    r  ^e 


2e       -|2K 


P   dt    .(3-77) 


Performing  the   integration,   we  have 


/ 


[-2y(t)-yM-yL 


yM   yL 


2K 


dt     < 


2e 
1 SL 


yM  yL- 


2K 


(T-to)    .(3.78) 


In  view  of  inequality  (3«69),  and  the  fact  that  yM  >  yT , 
it  follows  that 


Limit 

K  ->-°o 
P 


2e 


1  - 


yM  yL- 


2K 


P  . 


=  0  . 


(3-79) 


Observing  inequality  (3-78),  we  obtain 


Limit 

K  -►«> 
P 


/[ 


'rr2y(t)-yM-yLn 


yM  yL 


2K 


dt   =   0   . 


(3-80) 


The  theorem  is  proved. 
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IV.   THE  ALGORITHM 

This  section  describes  the  algorithm  used  for  minimizing 
the  augmented  performance  measure. 

A.   GENERAL  MINIMIZATION  STRATEGY 

1.  Sequence  of  Unconstrained  Sub-problems 

A  sequence  of  unconstrained  sub-problems  is  solved 
by  the  algorithm.   In  each  sub-problem  the  algorithm  mini- 
mizes the  augmented  performance  measure  for  given  values 
of  the  weighting  factors  (e  and  r)  and  the  inequality  con- 
straint penalty  term  power  (K  ) .   After  an  appropriate 
stopping  criterion  is  satisfied,  e  is  reduced,  K  is 
increased,  and  a  new  sub-problem  is  commenced  using  the 
optimal  solution  to  the  last  sub-problem  as  a  first  guess. 
This  procedure  is  repeated  until  enough  sub-problems  are 
completed  to  meet  a  second  stopping  criterion. 

The  algorithm  is  programmed  to  do  one  sub-problem 
on  each  computer  run.   The  results  are  stored  on  an  external 
storage  device  between  computer  runs  and  are  retrieved  at 
the  commencement  of  the  next  run  (new  sub-problem) . 

2.  Minimization  Strategy 

The  algorithm  minimizes  by  either  the  FNR  or  MNR 
method.   The  user  must  decide  which  method  to  use  on  each 
iteration.   This  is  a  matter  of  experimentation,  especially 
for  the  first  two  or  three  sub-problems.   An  effective  pro- 
cedure is  to  run  the  sub-problem  once  using  the  MNR  equation 


54 


i   \ 


throughout  and  once  using  the  FNR  equation  throughout.   Prom 
these  results  an  effective  minimization  strategy  can  generally 
be  deduced  for  the  sub-problem.   Occasionally  further  exper- 
imentation is  required.   This  experimentation  points  out 
the  advantages  of  using  separate  computer  runs  for  each 
sub-problem.   Once  a  sub-problem  is  completed  and  the  results 
stored,  the  computation  does  not  have  to  be  redone  each 
time  an  experimental  run  is  made  in  the  next  sub-problem. 

B.   COMMENCING  THE  PROBLEM 
1.   Initial  Decisions 

Three  interrelated  decisions  must  be  made  to  begin 
a  problem.   First,  the  number  of  time  points  K  must  be  chosen. 
Second,  the  number  of  coefficients  M  for  each  state  and 
control  expansion  must  be  chosen.   The  same  number  of  coeffi- 
cients is  used  for  all  expansions  in  a  given  problem  in 
this  dissertation,  but  this  is  not  a  requirement.   From  a 
theoretical  standpoint  it  Is  desirable  to  use  a  large  number 
of  coefficients  and  time  points  to  insure  that  an  adequate 
approximation  of  the  optimal  control  and  state  trajectory 
is  obtained,  but  practically,  computer  time  and  storage 
requirements  limit  the  number  of  each.   The  computational 
penalty  for  using  a  large  number  of  coefficients  is  the 
more  severe  of  the  two  as  the  number  of  equations  in  (2.36) 
and  (2.37)  which  must  be  solved  is  equal  to  the  total  number 
of  coefficients  plus  the  number  of  free  end  conditions.   The 
solution  of  equation  (-2.36)  or  (2.37)  represents  a  considerable 
portion  of  the  overall  computer  time. 
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The  third  decision  Involves  the  initial  values  of 

e,  r,  and  K  .   The  weighting  factor  r  for  the  inequality 

constraint  penalty  terms  is  held  constant  throughout  the 

entire  problem.   A  satisfactory  value  used  in  all  problems 

in  this  dissertation  for  all  inequality  constraint  penalty 

terms  is  r  =  100.   An  initial  value  of  K  which  has  worked 

P 

well  in  all  problems   is  K     =4.      Larger  values   of  K     gener- 

P  p 

ally  cause  computer  overflow  in  the  first  sub-problem. 
With  these  values  chosen  there  exists  a  region  of  e's  for 
which  the  first  sub-problem  will  respond  to  an  appropriate 
minimization  strategy.   This  acceptable  range  of  starting 
e's  is  different  for  each  problem  but  is  in  the  range 


10~5  <  e  <  10"3 


for  all  problems  solved  herein.   Numerical  experimentation 
is  the  only  method  available  to  determine  an  acceptable 
starting  e.   There  is  no  theoretical  requirement  to  use 
the  same  value  of  e  for  each  state  equation  equality  con- 
straint term  in  the  augmented  performance  measure  or  the 
same  K  in  each  inequality  constraint  term,  but  the  use  of 

different  e's  and  K  's  has  never  been  required. 

P 

2.   Initial  Guess  for  the  Unknowns 

Once  the  above  three  initial  decisions  are  made,  an 
initial  guess  for  the  vector  of  unknowns  c  is  required.   The 
vector  c  includes  all  coefficients  and  free  end  conditions. 
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All  coefficients  are  set  equal  to  zero  initially  unless 
there  is  good  reason  to  make  a  different  choice. 

C.   ITERATION 

1.   Required  Vectors  and  Matrices 

The  states  and  controls  are  calculated  at  each  time 
increment  by  evaluating  the  functional  expansions.   The  w 
vector  defined  in  (2.15)  is  calculated  using  these  states 
and  controls.   Next,  the  gradient  matrix  (2.34),  the  aug- 
mented performance  measure  (2.16),  the  symmetric  matrix 

T  T 

CV  w)   (V  w)  .  and  the  vector  (v  w)   w  are  calculated. 
c~  c   c-  c  c~  c  ~c 

At  this  point  the  algorithm  begins  the  iteration 
process  with  either  the  MNR  or  the  FNR  method  depending  on 
the  value  of  a  flag  set  by  the  user  (the  method  selected 
is  based  on  the  iteration  number  being  performed).   If  the 
MNR  method  is  to  be  used,  equation  (2.37)  is  formed.   If 

the  FNR  method  is  called  for,  the  three-dimensional  array 

2    T 
(2.35)  is  calculated  and  the  symmetric  matrix  (V  w)   w 

C  ***  C  *w  c 

is  formed.   It  is  prohibitive  to  store  the  entire  three- 
dimensional  array,  but  a  feasible  alternative  is  to  multiply 
each  matrix  in  this  array  by  w  as  the  matrix  is  calculated 
and  store  the  resulting  column  vector.   Once  a  matrix  in 
the  three-dimensional  array  is  multiplied  by  w,  it  Is  no 
longer  required  by  the  algorithm.   The  next  matrix  in  the 
array  Is  calculated  and  stored  in  the  same  storage  locations 

used  by  the  first  matrix.   Only  the  symmetric  matrix 

?    T 
(V  w)   w  need  be  stored.   The  total  increase  In  storage 
c  ~  c  ~c 
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requirements  of  the  FNR  method  over  the  MNR  method  using 
this  computation  technique  is  less  than  10  percent  in  the 
problems  solved  herein.   It  is  also  imperative  in  terms  of 

computation  time  to  take  full  advantage  of  the  symmetry  of 

2    T 
the  matrix  (V  w)   w  .   Due  to  this  symmetry  it  is  necessary 
c  ***  c  ***  c 

to  calculate  only  one  column  of  the  first  matrix  in  the 
array,  two  columns  of  the  second  matrix,  and  n  columns  of 
the  n   matrix.   By  taking  advantage  of  the  symmetry  the 
average  time  for  each  FNR  iteration  is  approximately  twice 
the  time  for  each  MNR  iteration. 
2.   Solving  the  Linear  System 

At  this  point  equation  (2.36)  or  (2.37)  is  formed 
and  must  be  solved  for  Ac.   This  is  a  linear  system  of  the 
form 

A  x  =  b  (4.1) 

and  is  solved  in  the  algorithm  by  one  of  three  methods 
available  to  the  user.   They  are: 

a.  Gauss  elimination  with  improvement  by  residuals 
using  total  pivoting, 

b.  Gauss  elimination  with  improvment  by  residuals 
using  main  diagonal  pivoting  and  a  computation  technique 
which  capitalizes  on  the  symmetry  of  A,  and 

c.  Gauss-Seidel  iteration. 

In  the  problems  solved  herein  the  number  of  unknowns 
varied  from  37  to  7^.   In  spite  of  the  large  number  of 
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unknowns  involved,  the  elimination  methods  required  less 
computation  time  to  solve  the  linear  system  than  the  Gauss- 
Seidel  iteration  method.   It  was  observed  for  the  problems 
solved  that  total  pivoting  was  not  required  In  the  elimina- 
tion method.   Method  b,  therefore,  was  the  most  economical 
and  effective  method  for  solving  the  linear  system  and  was 
used  in  all  problems.   Method  c  is  retained  in  the  event 
that  the  algorithm  is  used  to  solve  problems  with  a  larger 
number  of  unknowns. 

In  each  solution  of  equation  (4.1)  one  improvement 
is  made  using  residuals.  That  is,  after  equation  (4.1)  is 
solved, 

A  x  -  b  =  r  (4.2) 

is  formed.   The  system 

A  y  =  r  (4.3) 

is  solved  and  the  resulting  y  is  subtracted  from  x  to  form 
the  final  solution  to  equation  (4.1). 
3.   Interpolation 

Tabular  functions  of  two  independent  variables  are 
used  extensively  in  the  problems  to  represent  aircraft  and 
missile  parameters  accurately.   Parabolic  interpolation  is 
used  to  obtain  the  functional  values  in  these  tables  and 
the  required  first  and  second  partial  derivatives.   The 
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derivation  of  the  necessary  difference  equations  for  para- 
bolic interpolation  in  two  independent  variables  is  presented 
in  Appendix  C. 

Excerpts  from  the  tabular  data  used  in  the  problems 
is  presented  in  Appendix  B  along  with  graphical  representa- 
tions of  the  data.   The  data  represents  typical  supersonic 
aircraft  and  missile  performance  parameters  and  has  been 
obtained  from  several  sources.   Considerable  effort  was 
expended  to  smooth  the  data  before  the  tables  were  constructed 
since  finite  difference  methods  were  used  not  only  for 
functional  values  but  also  for  first  and  second  partial 
derivatives . 

4.   Stopping  Criteria 

Once  the  linear  system  is  solved,  a  new  c  vector  is 
calculated  from 


ci+1  =  c1  +  Ac1  .  (4.4) 


At  this  point  a  stopping  criterion  is  tested.   If 

|J  i  -  J0i+1|  :  ST0P1  ,  (4.5) 

1  a     a     — 

the  sub-problem  is  finished.   Otherwise  the  iteration  process 
is  continued.   At  the  completion  of  the  sub-problem  the 
results  are  stored  off  line.   The  computer  run  is  complete. 

To  begin  a  new  sub-problem  a  new  computer  run  is 
initiated,  recalling  the  results  stored  from  the  last 
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sub-problem.   Epsilon  is  decreased.  K  is  increased,  and 

P  * 

the  minimization  strategy  is  altered  by  the  user  as  required. 
Typically,  e  is  divided  by  a  factor  of  between  two  and  ten 
and  K  is  increased  by  two  or  four.   That  is, 


K  i+1  =  K  1   + 


2 

or 

4 


(4.6) 


More  ambitious  policies  usually  result  in  failure  of  the 
algorithm. 

Sub-problems  are  solved  until  a  second  stopping 
criterion  is  satisfied.   Several  criteria  are  possible  to 
end  the  problem.   A  method  used  successfully  involves 
observing 


Js*  +  Jp*  (iK7) 


and  stopping  when  this  sum,  which  represents  the  penalty 
terms  due  to  the  equality  and  inequality  constraints  without 
weighting  factors,  ceases  to  decrease  significantly  between 
sub-problems . 

5.  Flow  Chart 

A  flow  chart  of  the  algorithm  is  given  in  Figure  5 . 

6 .  Integration 

At  the  completion  of  the  last  sub-problem  a  check 
on  the  degree. of  satisfaction  of  the  state  equations  is 
obtained  by  comparing  the  state  expansions  with  the  state 
trajectory  obtained  by  integrating  the  state  equations  with 
the  control  expansions. 
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Read  constants,  flags,   initial   guess  for  c,  initial   conditions, 
weighting  factors,   inequality  constraint  power, 
fixed  final   conditions,  stopping  criteria. 


Evaluate  functional  expansions  at  each  time  point 


Calculate  w,  Vcw 


Form  ww,   (v  w)  w,  and  A  =  (v  w)   (v  w) 


Calculate  V    w 


i 

MNR 

?     T 
Form  (v"^  w)  w,  and 

A  =  A  +  (vc2w)Tw 

\ 

1 

Solve  for  Ac  = 

-  A_1(Vrw)Tw 

Calculate  new  c  vector,  c1       =  c    +  Ac 


Sub-problem  done.     Reduce  e's  and  increase  K  's 


No 


Figure  5 
Algorithm  Flow  Chart 
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V.   A  MISSILE  INTERCEPT  PROBLEM 

In  this  section  a  short-range  missile  intercept  problem 
is  solved.   An  air-to-air  missile  launched  from  an  attacking 
airplane  is  to  intercept  a  constant-velocity  target  in 
minimum  time.   The  missile  is  restricted  to  move  in  a  plane. 
The  orientation  of  this  plane  is  defined  in  three- 
dimensional  space  as  the  plane  containing  the  position  of 
the  missile  at  launch,  the  position  of  the  target  at  launch, 
and  the  velocity  vector  of  the  target.   The  assumptions 
applied  to  the  problem,  the  coorinate  systems  used,  the 
nomenclature,  and  the  derivation  of  the  equations  of  motion 
are  presented  in  Appendix  A. 

A.   PROBLEM  FORMULATION 
1.   State  Equations 

The  state  equations  derived  in  Appendix  A  are 


gTh  gpSa  2         gpSa  2         gW 

M  =  -rj-±±  Th  cosa  -  2w—  M  C  cosa  -  ^—   M  C^sina  -  -^  sine 
e  e  e  e 

(5.D 

gThM  Th  sina   gpSa  gpSa  gW   cos8 

6  =  wf  -is 2vr-  MCAslna  +  2TF-  MCNcosa  -  w;  — 
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(5.2) 


X  =  aMcose  -  aMTcosY  (5.3) 


Y  =  aMsine  -  aMTsiny    .  (5.*0 


The  states  are  Mach  number  M,  flight  path  angle  6,  relative 
range  X,  and  relative  cross-range  Y  as  defined  in  Figure  40. 
The  control  is  angle  of  attack  a. 

The  normalized  missile  thrust  Th  is  considered 
constant  until  missile  engine  burnout  tR  and  zero  thereafter. 


That  is, 


Th  =  1      ,  t  e  [0,tB]  (5.5) 


Th  =  0      ,  t  e  (tB,T]   .  (5-6) 


Other  parameters  considered  constant  are 


g 

32.1725 

ft. /sec.2 

Th   = 

3500 

lbs. 

W 
e 

200 

lbs. 

S 

0.35 

ft.2 

*B   = 

7 

sec . 

a   = 

1077. 

8 

ft. /sec. 

(5.7) 


p   =  0.001756   slugs/ft.3 


The  values  of  the  constants  in  equations  (5-7)  are  based  on 
an  engagement  at  10,000  feet  altitude.  It  is  convenient  to 
group  these  constants  as 


gTh 


C,  ■ 


1   aW, 


(5.8) 
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gpSa 
C2  =  2W—  (5.9) 


e 


C3  =  J-  •  (5.10) 


In  addition,  the  computer  solution  is  aided  by  normalizing 
the  states  X  and  Y  by  defining 

A  X 
x  =  £  (5.11) 


A  Y 
y  =  y  (5.12) 


where  X  and  y   are  assigned  nominal  values  of  10,000  feet 
With  these  adjustments  the  state  equations  are 


p  p 

M  =   C-jThcosct  -   C2M   C.cosa  -   C~2IVrCNsina   -   C  W   sine  (5-13) 


T.  ,  w  cose 

6  =  C±    in  ^lna  -  C2MCAsina  +  C2MCNcosa  -  0^   -^ (5-1^) 


x  =  ^  M  cose  -  §  Mtcosy  (5.15) 


y  =  y  M  sine  -  j   MTsinY    .  (5-16) 

2.   Tabular  Functions 

The  axial  and  normal  force  coefficients  C.  and  CM 
are  given  in  Appendix  B  as  tabular  functions  of  Mach  number 
and  angle  of  attack.   This  data  is  based  on  a  typical  air-to- 
air  missile. 
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3.  Performance  Measure 

The  performance  measure  for  this  problem  is 

J  =   /  dt   .  (5.17) 

(T 

4.  Inequality  Constraints 

Pour  inequality  constraints  are  required.   The  angle 
of  attack  must  satisfy 


-  f  <  a  <  f  .  (5.18) 


From  structural  considerations  the  load  factor  must  satisfy 


-50  <  |  (6M)  <  50   .  (5.19): 


5.   End  Conditions 

In  order  to  describe  the  initial  conditions  for  the 
problem  it  is  necessary  first  to  pose  the  problem  in  the 
(X,Y,Z_)  coordinate  system  shown  in  Figures  4l  and  42  of 
Appendix  A.   The  problem  chosen  for  presentation  in  this 
section  involves  a  target  in  a  shallow  climb  at  short  range 
with  a  slight  altitude  disadvantage  crossing  the  attacker's 
flight  path  extension  at  90°;  i.e. 

J^  =  15000  ft.                 (5-20) 

hT  =  -3000  ft.                (5.21) 

BT  =  90°  (5-22) 

6T  =  10°  (5.23) 
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Following  the  procedure  outlined  in  Appendix  A,  the 
remaining  unknown  parameters  are 

Hp        =   0.8  (5.24) 

Y  =   42.35°  (5.25) 

VT  =   51.526   lbs.  (5.26) 

C 

The  initial  conditions  as  computed  by  this  procedure  are 

M(0)  =  0.8  (5-27) 

6(0)  =  -49.573°  (5.28) 

x(0)  =  0  (5.29) 

y(0)  =  0   .  (5.30) 

The  final  conditions  as  computed  by  this  procedure  are 

x(T)  =  0.9920  (5.3D 

y(T)  =  -1.1645   .  (5.32) 

Note  that  x  and  y  are  the  normalized  relative  range  and 
relative  cross-range,  respectively,  as  defined  by  equations 
(5-11)  and  (5-12).   The  states  X  and  Y  in  equations  (5.11) 
and  (5.12)  are  defined  as  the  position  of  the  missile  in 
the  (X,Y)  coordinate  system  shown  in  Figure  40.   At  t  =  0 
the  missile  is  at  the  origin  of  the  (X,Y)  system,  hence 
equations  (5.29)  and  (5-30)  apply.   At  t  =  T  the  missile 
must  intercept  the  target.   Since  the  (X,Y)  system  becomes 
fixed  with  respect  to  the  target  at  launch,  x(T)  and  y(T) 
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are  equal  to  the  normalized  coordinates  of  the  target  in  the 
(X,Y)  system  at  launch  and  are  given  by  equations  (A. 93)  and 
(A. 95). 

B.   THE  EPSILON  METHOD  FORMULATION 

1.   The  Augmented  Performance  Measure 

Using  the  penalty  functionals  described  in  Section 
III  for  inequality  constraints,  the  augmented  performance 
measure  is 

J*=   -fdt 

T 

+  ly      m  _  Ojlhcosa  +  C2M2CAcosa  +  CgM^sina  +  C-W  sine]     dt 
T 

i  C  r  •       th  i  w  cose 1 2 

+  eJ     L6  "Cl'^2L  +  C2MCASijla  "  C2MCNC0Sa  +  C3  ~V~J    dt 


0 

T  2 

+  j J     \x  -  j  Mcose  +  |  Mpcosy]     dt 

/  T  2 

+  I  J    [y  ~  7  ^ih8  +  7  V1^]   dt 


(5.33) 


rT   m  n2K         /-  Tr  AM  -]2K 

♦'.(  ft]  p*=  +  <  [§&]  p* 


where  e  and  r  are  weighting  factors  and  6  is  a  constant  used 
to  make  minor  adjustments  in  the  admissible  regions.   In  all 
problems  6  is  given  a  value  of  1.03.   This  adjustment  is 
applied  to  all  admissible  regions  of  constrained  states  and 
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controls.   The  power  K  is  limited  in  computation  to  approxi- 
mately 30  before  computer  exponent  overflow  problems  develop. 

With  this  value  of  K  it  is  desirable  to  adjust  all  admis- 

P 

sible  regions  slightly  as  can  be  seen  from  Figure  6. 


r2y-yM-yLi60 


Figure  6 
Adjusted  admissible  regions 
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The  required  elements  of  w  are 


w 


k   =    [«k   "    ClThcosak  +    C2\2\coa«*   +    C2Mk2Vina* 


+   C3WcBinek]  [&\k        ,        k  =   1,2,..., 


(5.34) 


K 


wK+k =  [K  -  °1 


Thsina, 


T£~   +  C2MkCAksinak  "  C2MkCNkCOSak 


(5.35) 


w  cose, 

c    k 


"c  — k]  \LtlH  k  =  1  2 


K 


w 


2K+k 


-  [xk  -  |  Mkcosek  +  |  MTcosY]  [^]%  ,  k=l,2,...,K   (5.36) 


w 


3K+k 


=  [*k  -  7  Mksln9k  +  7   MTsinY]  [t"]%  '  k=1»2»-.-»K   (5.37) 


w 


4K+k  = 


20^1 
tt6 


K 


P  (rAt)% 


IV   ""    1  j  £.  j  «  «  i  j  1\ 


(5.38) 


w 


5K+k 


raekMk 
50g6 


K 


p(rAt)%    ,   k  =  1,2,.. ,,K 


(5.39) 


W6K+1  =  C'CK-l)At] 


(5. HO) 


where 


ML  =  M(t) 
k 


t=(k-l)At 


etc. 
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2.      Functional  Expansions 

The   state   and   control   expansions  written  In  discrete 
form  are 


\  =  M1  +  JS_±  (k-1)    +    £    amsin  mTT^i        »    k-l,2,...,K      (5. Hi) 

m=l 

e„-e,  m  fir  i\ 

6k  =    6l   +  -T^r  (k_1)    +    £     Vin        K-i        »    k=l,2,...,K      (5.12) 

m=l 

Xk  =   Xl   +  -fcT1   (k"1)    +     £     cmSln  m7TK-l1)    "    k=l,2,...,K      (5.43 

m=l 

•  y^-yn  M        imr(k-l)   .  ,  9     „   ,R  y, 

yk  =  yl  +  "fa1  (k"1)  +  E  dmsin   K_1    '  k-1»2>-"'K   C5-Z,4) 

m=l 

ak  =  al  +  "TEXT   (k_1)  +  ^  emsln  ""'K-I1''  ■  k=l,2,...,K.  (5.^5) 

m=  1 


3.   Vector  of  Unknowns 

The  elements  of  the  vector  c  are 


T 
c   =  w^ »&2»  •  *  •  *a]yi»  i»°2 ' ' ' ' '  M,cl* c2* '  *  *  ,CMJ 


d1,d2, . . . ,dM,e1,e2, . . . ,eM,MK,eK,a1,aK,At)  . 


(5.^6) 


1.   Partial  Derivatives  of  the  Tabular  Function 

The  values  of  C.  and  C.T  are  obtained  from  the  tables 

30A    =>CA 
by  parabolic  interpolation  along  with  the  values  of  -ttt—   ,  — —  , 

a  vl  3  ct 
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>\ 

92CA 

9CN 

9CN 

>\ 

*\ 

*  2     » 

'   » 

* 

i 

2  '    » 

2       ' 

3cT 

3M3a 

3M 

8a 

3M 

3ci 

32CA  __ 

5   >    »!   »   — - ~~ ~~  *      i   ~"~" ™  t    o   >    i   i   and 
9M 

32CN 

A„a   .    The  procedure  is  outlined  in  Appendix  C. 

5 .   First  Partial  Derivatives 

The  first  partial  derivatives  indicated  in  equation 
(2.3*1)  are  required.   These  partial  derivatives  are  easily 
obtained  from  equations  (5.3*0  thru  (5.*J0)  by  taking  the 
partials  of  these  expressions  with  respect  to  the  vector  c. 
The  expressions  are  too  numerous  to  include.   A  typical  term 
is 


^   3wk  3M*   3W£  3\ 

9am~  ^k"9am   9\  9am   '  ^   ^ 


For  1  <  k  <  K,  w,  is  given  by  equation  (5.3*0  and  the  partial 
derivatives  indicated  above  are 


3w 

3M 


k  =  |"Atl%  (5i48) 

i.    L  e  J 


3*,    r  o  9CA 


3MJ  =  [2C2MkCAkCOSak  +  C2Mk2  3M^  C0Sak  +  2C2MkCNkslnak 

(5.*»9) 
2wk      3M. 


9CN 
2        1Nk 


-    '     —  ^][ff 


k 


3M 

^-TK^m^^T1  (5.50) 
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k  _  o-»n  nnr(k-l) 

9a  '  sln   K-l 
m 


(5.51) 


Notice  that  C.   and  CM  are  functions  of  M.  as  well  as  a,  . 
Ak      \  K  K 

6.   Second  Partial  Derivatives 

The  second  partial  derivatives  indicated  in  the  three- 
dimensional  array  (2.35)  are  also  required.   Again  the  expres- 
sions are  too  numerous  to  list.   A  typical  term  for  1  <  k  <  K 
is 


2 
3  w 


k   _  3 


raw, 


9Vam   9a* 


8a_ 


m  -1 


(5.52) 


3W  A 
Letting  — ^  =  R 
9am 


we  have 


3   w 


k        _   3R 


9R      8Mk   +    3R_  *\ 


9Vam       9a*        9%  8a£        3Mk  3a£      ' 


(5.53) 


The  partial  derivatives  indicated  above  are 


3R 

3^k 


=  0 


(5.5*0 


3R_ 

3M, 


8CA 

2C2CAkCOSak  +  *C2Mk  3Mk 


2 

9  CA 

k  2     k 

—  cosa,  +  C0M.   ~—  cosol 

K    d   K   3M,        K 


+  2C~C.T  sina,  +  4C0M 

2  Nk    k     2  .  _k 


2 

9CN  9  CN 

-k  wr  slnak +  c2\  rnr 


sina 


3M, 


k 


[f]! 


^sin  mTr(k-l) 


K-l 


(5.55) 
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8Mk         &TT  gTT(k-l)  fr-     cC\ 

9i7  =  (K-l)At  COS    K-l  (5*56> 


^=  sin  M^i   .  (5.57) 


C.   RESULTS 

The  problem  was  solved  twice:   once  using  8  coefficients 
for  each  expansion  (problem  A)  and  once  using  12  coefficients 
for  each  expansion  (problem  B),   In  both  cases  K  =  21  time 
points  (20  time  intervals)  were  used.   The  initial  guess 
for  the  c  vector  was:  i 


all  expansion  coefficients  =  0 
M(T)  =  1.4 

a(0)  =  20°  15O0J 

a(T)  =  0° 

At  =  7/20  sec.  (T  =  7  sec.) 


1.   Problem  A  -  8  Coefficients  for  each  Expansion 

Pour  sub-problems  were  required  to  solve  the  problem. 
Table  3  gives  the  weighting  factor  values,  optimization 
strategy,  and  computer  time  for  each  sub-problem.   Table  4 
gives  the  components  of  the  minimum  augmented  performance 
measure  for  each  sub-problem  where 


Ja  *  J  +  Fs  +  rJP  (5'59) 


7^ 


\      \ 


sub- 
problem 

e 

R 

K 
P 

I* 

optimization 
strategy** 

C.P.U. 
time 

1 

1.0  x  10"5 

100 

4 

8 

MMMMMMMM 

2'28" 

2 

0.67  x  10"5 

100 

6 

6 

FFFFFF 

3' 38" 

3 

0.5  x  10"5 

100 

8 

2 

FF 

1*54" 

k 

0.5  x  10"5 

100 

32 

1 

F 

1'32" 

*  number  of  iterations  required 
**  M  -  MNR  method 
F  -  FNR  method 


Table  3 


Weighting  factors,  optimization  strategy,  and  C.P.U.  time  for 
missile  intercept  problem  A. 


sub- 
problem 

V 

J* 

* 

J 
s 

* 

JP 

1 

16.95 

7.538 

0.6853  x  lO"4 

0.2554  x  10"1 

2 

14.52 

7.585 

0.4626  x  10"1* 

0.1287  x  10~9 

3 

16.83 

7.585 

0.4624  x  10"4 

0.4619  x  10"13 

4 

I6.83 

7.585 

0.4624  x  10 

0.3691  x  10"53 

Table  4 

Components  of  the  minimum  augmented  performance  measure  for 
missile  intercept  problem  A. 
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Figure  7  is  a  plot  of  the  augmented  performance 
measure  vs.  iteration  number. 
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Figure  7 

Augmented  performance  measure  vs.  iteration 
number  for  missile  intercept  problem  A 
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Iterations  performed  by  the  PNR  method  are  indicated  by  a 
solid  line.   Iterations  performed  by  the  MNR  equations  are 
indicated  by  a  broken  line.   The  figure  shows  several  signi- 
ficant characteristics: 

a.  the  failure  of  the  FNR  method  on  the  first  itera- 
tion of  the  problem  (the  initial  guess  is  too  far  from 
optimum  to  allow  the  FNR  method  to  converge); 

b.  the  superior  terminal  convergence  produced  by 
the  FNR  method  close  to  the  minimum; 

c.  the  oscillatory  results  produced  by  the  MNR 

method  as  the  minimum  is  approached.   After  the  commencement 

of  sub-problem  2  as  shown  in  Figure  7,  the  FNR  method  is 

used  exclusively  to  the  end  of  the  problem.   For  comparison, 

sub-problem  2  is  commenced  with  the  MNR  method  and  allowed 

to  run  to  23  iterations.   At  the  beginning  of  sub-problem  3» 

J  increases  slightly  and  never  returns  to  the  minimum  value 

obtained  in  sub-problem  2.   This  is  an  indication  that  e  has 

reached  a  point  where  smaller  values  have  little  influence 

on  the  reduction  of  J  . 

s 

Figure  8  is  a  plot  of  the  performance  measure  vs. 

iteration  number.   The  performance  measure  (final  time) 

increases  with  iteration  number.   As  the  algorithm  iterates, 

At  is  being  increased  to  reduce  J   (the  term  in  the  augmented 

s 

performance  measure  reflecting  the  degree  of  non-satisfaction 
of  the  state  equations)  while  insuring  that  constrained 
states  and  controls  are  kept  within  admissible  bounds  by 
reducing  or  holding  down  J  .   With  small  values  of  e  and 
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Sub-problem  number 
1 


3  4 


Figure  8 

Performance  measure  vs.  iteration  number 
for  missile  intercept  problem  A. 


large  values  of  K  the  end  result  of  minimizing  the  augmented 
performance  measure  is  a  control  history  and  state  trajectory 
which  minimizes  the  time  to  intercept  while  satisfying  the 
state  equations  and  inequality  constraints;  all  to  a 
reasonable  degree  of  accuracy. 
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Figure  9  is  a  plot  of  J  and  J  vs.  iteration 

s  p 


number. 
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Sub-problem  number 
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16    18 


Figure  9 

J  and  J  vs.  iteration  number 
s      p 

for  missile  intercept  problem  A. 


After  K  is  increased  at  the  beginning  of  sub-problem  2,  the 
inequality  constraint  penalty  terms  J  become  very  small. 
This  is  because  the  angle  of  attack  is  well  within  its 
admissible  region. 


79 


It  is  evident  from  Figures  7>  8,  and  9  that  two 
sub-problems  are  sufficient  to  obtain  a  reasonable  solution. 
The  overal  stopping  criterion 


(Js*+JpV  -   (Js*+JpV+1|  <10"6  (5.60) 


where  I  is  the  sub-problem  number  is  satisfied  after  the 

third  sub-problem.   However,  at  this  point  the  value  of  K 

p 

is  8  which  is  not  large  enough  to  provide  desirable 
penalty  functionals  for  the  inequality  constraints.   It  is 
necessary  to  provide  as  large  a  value  of  K  as  the  computer 
will  allow  for  the  final  sub-problem  to  insure  that  the 
minimization  is  not  influenced  by  the  inequality  constraint 
penalty  terms  when  the  constrained  states  and  controls  are 
within  their  admissible  regions.   Accordingly,  a  final  sub- 
problem  is  performed  with  K  =  32.   The  algorithm  is  able 
to  handle  the  increase  in  K  from  8  to  32  in  one  step  only 
because  no  constraint  boundaries  are  active  in  the  solution. 

Figure  10  is  a  plot  of  the  angle  of  attack  expansion 
computed  at  the  end  of  the  last  sub-problem.   The  region  of 
admissible  angles  of  attack  is  shown. 

Figures  11  and  12  are  plots  of  Mach  number  and 
flight  path  angle  vs.  time.   In  each  plot  two  curves  are 
shown;  one  is  the  expansion  for  the  state  as  computed  at  the 
end  of  the  last  sub-problem;  the  other  is  the  state  trajec- 
tory obtained  by  numerically  integrating  the  state  equations 
with  the  optimal  control  expansion. 
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Figure  10 

Angle  of  attack  vs.  time  for 
missile  intercept  problem  A. 
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Figure    11 

Mach  number  vs.    time    for 
missile   intercept   problem  A. 
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Figure  12 

Flight  path  angle  vs.  time  for 
missile  intercept  problem  A. 
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Figure  13  is  a  plot  of  relative  range  X  vs.  relative 
cross-range  Y.   As  before  two  curves  are  presented:   one 
represents  the  expansions  of  the  states  as  computed  at  the 
end  of  the  last  sub-problem;  the  other  is  the  state  trajec- 
tory obtained  by  numerically  integrating  the  state  equations 
with  the  optimal  control  expansion. 


Relative  range  (ft.  x  10  ) 


10 


12 


integration 


Figure  13 

Relative  range  vs.  relative  cross-range 
for  missile  intercept  problem  A. 
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Figure  14  is  a  plot  of  range  X'  vs.  cross-range  Y1 
where  both  quantities  are  obtained  by  transforming  the 
expansions  of  the  states  X  and  Y  obtained  at  the  end  of  the 
last  sub-problem  from  the  (X,Y)  coordinate  system  to  the 
inertial  coordinate  system  (X',Y')  fixed  at  the  missile 
launch  point  (Figure  *J0). 
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Figure  14 

Range  vs.  cross-range  for 
missile  intercept  problem  A. 
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Figure  15  is  a  plot  of  load  factor  vs.  time  where 
the  load  factor  is  given  by 


n  =  -  (9  M) 
g 


(5.61) 


and  the  states  used  in  equation  (5.6l)  are  the  state  expan- 
sions obtained  at  the  end  of  the  last  sub-problem. 


Figure  15 

Load  factor  vs.  time  for 
missile  intercept  problem  A. 
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Although  load  factor  Is  not  a  state  but  a  function  of  states, 
the  plot  is  important  because  it  shows  that  the  load  factor 
constraints  as  given  by  (5.19)  are  not  active. 

It  might  be  suspected  that  the  optimal  trajectory 
is  a  maximum  performance  turn  limited  by  a  constraint 
boundary  (angle  of  attack  or  load  factor)  followed  by  a 
straight  line  path  to  intercept.   This  is  not  the  case.   The 
initial  turn  rate  of  the  missile  is  small  compared  to  its 
maximum  turn  rate  capability.   This  is  due  to  the  high 
induced  drag  associated  with  high  angle  of  attack  turns 
which  would  reduce  the  missile's  longitudinal  acceleration 
capability.   Also  there  is  no  straight  line  segment  to  the 
trajectory  although  the  turn  rate  of  the  missile  is  very 
small  at  intercept. 

The  first  sub-problem  was  also  solved  with  an 
initial  guess  for  the  c  vector  of: 


all  expansion  coefficients  =  0 
M(T)  =  2.5 

6(T)  =  0°  (5.62) 

a(0)  =  10° 
a(T)  =  0° 

At  =  8/20  sec.  (T  =  8  sec.) 


The  sub-problem  reached  a  minimum  of  7.5^4  which  compared 
favorably  with  the  minimum  reached  by  the  first  initial 
guess  given  by  equations  (5.58).   This  gives  an  indication 
that  the  minimum  attained  is  the  global  minimum. 
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2.   Problem  B  -  12  Coefficients  for  each  Expansion 
Four  sub-problems  were  required  to  solve  this 
problem  also.   Tables  5  and  6  summarize  the  performance  of 
the  algorithm. 


sub- 
problem 

e 

r 

K 
P 

I* 

optimization 
strategy** 

C.P.U. 
time 

1 
2 

3 
4 

1.0  x  10"5 
0.67  x  10"5 
0.5  x  10~5 
0.5  x  10~5 

100 
100 
100 
100 

4 

6 

8 

32 

8 
2 
2 
2 

MMMMMMMM 

PP 

FF 

FF 

3'  34" 
2'20" 

2'19" 
2120" 

*  number  of  iterations  required 
**  M  -  MNR  method 
F  -  FNR  method 

Table  5 

Weighting  factors,  optimization  strategy,  and 
C.P.U.  time  for  missile  intercept  problem  B. 


sub- 
problem 

V 

J* 

J  * 
s 

J  * 

P 

1 

11.51 

7.537 

0.2125  x  10~4 

0.1849  x  10"1 

2 

9.69 

7.612 

0.1385  x  10"2* 

0.3412  x  10"8 

3 

10.65 

7.628 

0.1513  x  10"4 

0.3474  x  10"11 

4 

10.31 

7.631 

0.1339  x  10_2j 

0.6139  x  lO"43 

Table  6 

Components  of  the  minimum  augmented  measure  for 
missile  intercept  problem  B. 
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Figure  16  is  a  plot  of  the  augmented  performance 
measure  vs.  iteration  number  and  corresponds  to  Figure  7 
for  problem  A. 
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Figure  16 

Augmented  performance  measure  vs.  iteration  number 
for  missile  intercept  problem  B. 
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A  second  failure  mode  of  the  MNR  method  close  to  the  minimum 
is  shown.   The  MNR  method  produces  a  divergent  J  in  the 
second  sub-problem. 

Figure  17  is  a  plot  of  range  X'  vs.  cross-range  Y' 
obtained  in  the  same  manner  as  in  problem  A  (Figure  14) . 


-8   3  -10 

Cross-range  (ft.  x  10  ) 


12 


Figure  17 

Range  vs.  cross-range  for 
missile  intercept  problem  B. 
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A  comparison  of  the  tables  and  figures  for  problems 


A  and  B  show  that  the  optimal  control  and  trajectory  have 
not  been  markedly  affected  by  the  increase  in  the  number  of 
coefficients  from  8  to  12  for  each  expansion.   It  is  prohib- 
itive in  terms  of  computation  time  and  storage  requirements 
to  increase  the  number  of  coefficients  further. 
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VI. .  A  CLIMB  PERFORMANCE  PROBLEM 

In  this  section  a  climb  performance  problem  is  solved. 
A  supersonic  fighter  aircraft  is  to  climb  from  sea  level  to 
high  altitude  in  minimum  time. 

Plight  test  experience  has  shown  that  to  climb  to 
altitudes  above  the  tropopause  in  minimum  time  a  supersonic 
fighter  must  execute  a  maneuver  which  typically  includes: 

a.  a  subsonic  climb  to  an  altitude  near  the 
tropopause; 

b.  a  level  or  near  level  acceleration  to  some 
supersonic  Mach  number; 

c.  a  "zoom"  climb  to  the  desired  altitude  trading 
kinetic  energy  for  potential  energy.   This  technique  has 
been  used  extensively  in  the  past  decade  for  establishing 
climb  records  and  in  fighter-interceptor  tactics  to  attain 
altitudes  higher  than  the  aircraft's  service  ceiling  for 
short  periods  of  time. 

Several  factors  contribute  to  the  optimality  of  this 
type  of  maneuver.   They  are: 

a.  a  fighter's  maximum  Mach  number  or  "placard" 
limit  which  arises  from  dynamic  pressure  and/or  thermal 
limitations  and  is  a  function  of  altitude; 

b.  a  fighter's  transonic  drag  characteristics; 

c.  air  temperature  variation  with  altitude; 

d.  air  density  variation  with  altitude; 

e.  turbojet  engine  maximum  thrust  variation  with 
altitude , 
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Optimization  techniques  were  first  applied  successfully 
to  this  problem  by  Bryson  [Refs.  1  and  2].   The  method  of 
steepest  descent  was  used  successfully  to  predict  the  type 
maneuver  described  above  for  a  typical  supersonic  aircraft. 

The  epsilon  method  is  applied  to  the  problem  herein  to 
demonstrate  the  method's  power,   A  direct  comparison  of 
methods  is  not  made  as  the  mathematical  model  used  here  has 
been  improved  considerably  over  that  used  in  Reference  2. 

A.  pr6blem  FORMULATION 
1.   State  Equations 

The.  state  equations  for  this  problem  derived  in 
Appendix  A  are 


M  =  £i£  cosa  -  |  siny g§ —  aM  CD  (6,1) 

e 

^  =   gTh   sina  +  f^aMCT    .   g£osy_  (g>2) 


M      2w~   w      L     aM 


h  =  ^   sin  y  .  (6.3) 

HL 


The  states  are  Mach  number  M,  vertical  flight  path  angle  y> 
and  normalized  altitude  h.   The  control  Is  angle  of  attack  a. 
Parameters  considered  constant  are  the  gravitational  constant 
g,  sea  level  standard  density  p  ,  aircraft  wing  area  S, 
aircraft  weight  W  ,  and  the  altitude  of  the  tropopause  under 
standard  atmospheric  conditions  H"L<   These  constants  are 
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g  =  32.1725   ft. /sec.2 

Po  =  0.002378  slugs/ft.3 

S  -  400  ft.2  (6.4) 

W  =  39,000  lbs. 

HL  =  36,089  ft. 

The  remaining  parameters  in  equations  (6.1)  thru 
(6.3)  vary  with  flight  and  atmospheric  conditions.   These 
variations  are  represented  in  either  tabular  form  or  by 
empirical  relations.   These  tabular  and  empirical  relations 
are  critical  to  the  problem  as  they  represent  the  mathemat- 
ical equivalents  of  the  factors  a  through  e  listed  in  the 
introduction  to  this  section. 

It  is  convenient  to  define  the  constant 


A  gp  S 

c"tt    •  (6'5) 

e 


With  the  definition  (6.5)  incorporated  the  state  equations 
are 


M  =  ^  cosa  -  §  siny  -  caaM2Cn  (6.6) 

a        a  D 


;  .  g|h  sina  +  oaoM0L  .  go^  (6.?) 


aM 


h  =  an  siny  t  (6§8) 
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2.   Empirical  Relations 

Empirical  relations  are  used  for  air  density  ratio  a, 
maximum  Mach  number  fi.,  and  the  speed  of  sound  a  as  functions 
of  normalized  altitude  h.   Air  density  ratio  and  normalized 
altitude  are  defined  by 


=  ■£-  (6.9) 

Po 


and 


h  =  g-  (6.10) 

HL 


These  empirical  relations  which  are  discussed  in  Appendix  D 
are  repeated  here  for  convenience.   They  are 


-c,h       -c9h 
a   =  e  ■   +  c3he  d  (6.11) 

MM  =  2.1  -  l.le"2*4h  (6-.12) 


a  =  aQ(l  -  c?h)   ,  h  <  1  (6.13) 


=  971  ft. /sec.  ,  h  >  1   .  (6.14) 


where 


c1  =  1.54100  (6.15) 

c2  =  1.80445  (6.16) 

c3  ■  0.4130  (6,17) 

c4  =  0.1331  .                          (6.18) 
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3.   Tabular  Functions 

In  situations  where  parameter  variations  cannot  be 
adequately  represented  by  empirical  formulas,  a  table  of 
values  is  used.   Tables  are  used  for  lift  and  drag  coeffi- 
cients as  functions  of  Mach  number  and  angle  of  attack  for 
a  typical  supersonic  fighter.   Excerpts  from  these  tables 
are  presented  in  Appendix  B. 

The  thrust  Th  appearing  in  equations  (6.6)  and  (6.7) 
is  normalized  maximum  thrust  as  it  is  assumed  that  since 
the  aircraft  must  climb  to  altitude  in  minimum  time,  its 
power  plant  will  always  be  operated  at  maximum  thrust. 
Maximum  thrust  is  normalized  with  respect  to  sea  level 
static  maximum  thrust  ThM  and  is  given  by 


Th 

Th  =  -JL  (6.19) 

1  M 
o 


where 


ThM  =  3^,000  lbs.  (6.20) 

o 


Maximum  thrust  is  given  as  a  tabular  function  of  Mach  number 
and  altitude  for  the  fighter  under  consideration. 
*J.   Performance  Measure 

The  performance  measure  for  this  problem  is 

T 

J  =  /   dt   .  (6.21) 

J0 
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5.   Inequality  Constraints 

The  following  state  and  control  inequality  constraints 
are  imposed: 


0  <  M  <  MM  (6.22) 


-6°  =  ctL  <  a  <  aM  =  24°   .  (6,23) 


The  maximum  Mach  number  M„  constraint  represents  the 
"placard"  limit.   M„  is  a  function  of  altitude  and  is  given 
by  an  empirical  relation  as  discussed  in  Section  VI. 2. 

The  angle  of  attack  constraint  aM  is  set  at  an 
angle  of  attack  slightly  above  that  for  maximum  lift  coeffi- 
cient.  The  minimum  angle  of  attack  a  is  set  slightly  below 
that  for  minimum  lift  coefficient  thus  simulating  aerodynamic 
stall. 

6.   End  Conditions 

The  initial  conditions  are 


M(0)  =  M  =  0.6 
o 


y(o)  =  0 


h(0)  =  hQ  =  0   . 


C6.24) 

(6.25) 

(6.26) 

The  final  condition  is 


h(T)  =  hp  -  60>°g°  ft>     .  (6.27) 


'L 
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B.   THE  EPSILON  METHOD  FORMULATION 

1.   The  Augmented  Performance  Measure 

Using  the  penalty  functional  described  in  Section  III 
for  inequality  constraints,  the  augmented  performance  measure 
is 


-f 


J*%     " 


+  -  f     [m  -  ■£^h-  cosa  +  ^  sinY  +  CaaM2CD]  dt 

♦  kf!  [i  -  if  — ]2  at 


(6.28) 


*.-jT«-P«  *  <[^P« 


where  e  and  r  are  weighting  factors,  6  is  a  constant  used  to 
make  minor  adjustments  to  the  admissible  regions,  and  d  is 
the  midpoint  of  the  admissible  angle  of  attack  region;  that 
is 


d  =  M  ,  L  (6.29) 
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The   required  elements   of  w  are 


wk'=[\- 


—~  cosa,    +  -&- 


■*"*  +  flV*%J  [f F  ■  (6,30) 


k=l,2, . . .,K 


w 


■[ 


gH^  sina^. 


gcosy. 


K=l j  2 , ... j  K 
a,  M, 


w 


2K+k 


(6.3D 


=    [\  ~  "T^  SlnYk]  [*T"]%      '      k=l,2,...,K  (6.32) 


f   2M,  IK 


w 


3K+k  °  [  MJJ 


-  1 


"p    (rAt)%      ,  k=l,2,..., 


K 


(6.33) 


w 


4K+k 


ak"d 
aM"d 


K 


p    (rAt)%      ,  k=l,2,...,K 


(6.34) 


W5K+1   =    t(K-l)At]% 


(6.35) 


2.      Functional  Expansions,    Unknowns,    and  Partial 
Derivatives 

The    state   and   control  expansions   are   of  the   same 

form  as   the  problem  in  Section  V  and  are  not   shown.      The 

elements   of  the    c  vector  are 


T 
c     =    (a,  ,a5|  •  •  .  fB.fj.fU-,  jDnj » •  •  >oM,  C-.  jC5,  •  ■ .  >cM, 


ll>a'2 ; 


*M»"1*"2 


M»wl»w2' 


'M; 


d1,d2>. . . ,dM,MK,YK,a1,aK,At) 
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(6.36) 


l       v 


where  the  a  's  represent  the  Mach  number  expansion  coeffi- 

cients,  the  b  ' s  represent  the  vertical  flight  path  angle 

coefficients,  the  c  's  represent  the  altitude  coefficients, 
'      m  * 

and  the  d  '  s  represent  the  angle  of  attack  coefficients, 
m 

The  first  and  second  derivatives  of  the  empirical 
relations  (6.11)  thru  (6.14)  are  required.   These  expressions 
are 


r\rt  -C-h  -Cph 

gg  =  -C;Le  x  +  (c3-c2c3h)e  (6.37) 

A2  ~  -c,h  -c«h 

2-2.  =  c/e  L     -   [(c_-cpc,h)c0  +  c?c„Je  d  (6.38) 

dh     x  5  J  d         d   J 

^=2.64e-2'4h  (6.39) 

i^  =  -6.336e-2<I,h  ■     (6.40) 

86-  "aoc7   »   h  <  1   ..  (6.41) 


=   0   .    ,   h  >-l  (6.42) 


A  _ 
dh 


,-  =  0   .  (6.43) 


The  first  and  second  partials  of  the  tabular  functions 
for  lift  coefficient,  drag  coefficient,  and  maximum  thrust 
with  respect  to  their  independent  variables  and  the  elements 
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of  w  given  by  equations  (6. 30)  thru  (6.35)  with  respect  to 
c  are  obtained  in  the  same  manner  as  in  the  problem  in 
Section  V. 

C.   RESULTS 

Two  problems  were  solved.   In  problem  A  the  aircraft  was 
to  climb  from  sea  level  to  60,000  feet  in  minimum  time. 
Problem  B  encompassed  a  series  of  problems.   The  results  of 
problem  A  were  used  as  a  first  guess  for  the  solution  of  a 
minimum-time-to-climb  profile  from  sea  level  to  61,000  feet, 
which  in  turn  was  used  as  a  first  guess  for  a  climb  to 
62,000  feet,  etc.   In  this  manner  optimal  control  and  state 
trajectories  were  obtained  for  minimum-time-to-climb  profiles 
from  sea  level  to  altitudes  from  60,000  to  70,000  feet  in 
thousand-foot  increments.   In  both  problems  8  coefficients 
for  each  expansion  and  4l  time  points  ( 40  time  intervals) 
were  used. 

1.   Problem  A  -  A  Climb  from  0  to  60,000  Feet 
The  initial  guess  for  the  c  vector  was: 

all  expansion  coefficients  =  0 
M(T)  =  0.9 

Y(T)  =  45°  ,g  nil) 

a(0)  =  1°  tt>.44; 

a(T)  =  1° 

At  =  120/40  sec.  (T  =  2  min.) 

Four  sub-problems  were  required  to  solve  the  problem. 
Tables  7  and  8  summarize  the  performance  of  the  algorithm. 
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sub- 
problem 

e 

r 

K 
P 

I* 

optimization 
strategy** 

C.P.U. 

1 

0.2  x  10_/| 

100 

4 

9 

FMMMFFFFF 

3' 33" 

2 

-4 

0.1  xlO 

100 

6 

9 

FMMMMMMFF 

3*24" 

3 

0.67  x  10"5 

100 

8 

8 

MFFFFFFF 

3' 38" 

4 

0.5  x  10"5 

100 

32 

3 

FFF 

2'59" 

*   number  of  iterations  required 
**   M  -  MNR  method 
F  -  FNR  method 


Table  7 


Weighting  factors,  optimization  strategy,  and  C.P.U. 
time  for  climb  performance  problem  A. 


sub- 
problem 

Ja* 

J* 

V 

V 

1 

422.  5 

129.1 

0.4496  x  10~2 

0.6855 

2 

442.2 

199.3 

0.2154  x  10"2 

0.2742 

3 

489.3 

238.6 

0.1544  x  10~2 

0.1904 

4 

451.5 

258.0 

0.9677  x  10"3 

O.8838  x  10' 

-12 

Table  8 

Components  of  the  minimum  augmented  performance  measure 
for  climb  performance  problem  A. 
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Neither  the  angle  of  attack  or  maximum  Mach  number 
constraints  are  active  in  this  problem.   The  largest  Mach 
number  attained  in  the  climb  is  1.53  at  an  altitude  of 
23,000  feet.   Consulting  Appendix  D,  the  maximum  Mach  number 
at  this  altitude  is  1.88.   However,  since  both  constrained 
parameters  approach  their  boundaries,  a  large  value  of  K 
is  required  to  obtain  small  J  contributions  to  the  aug- 
mented performance  measure. 

Figure  18  is  a  plot  of  the  augmented  performance 

measure  vs.  iteration  number. 
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Figure    18 

Augmented  performance  measure   vs.    iteration 
number   for   climb   performance   problem  A. 
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Iterations  performed  by  the  FNR  method  are  indicated  by  a 
solid  line.   Iterations  performed  by  the  MNR  method  are 
indicated  by  a  broken  line.   As  observed  in  the  previous 
problem,  the  FNR  method  results  in  excellent  terminal 
convergence .   The  MNR  method  performs  well  when  the  c 
vector  is  far  from  optimum  as  indicated  by  relatively  large 
augmented  performance  measures.   At  the  commencement  of 
sub-problem  1  the  MNR  method  fails  presumably  because  the 
initial  guess  for  the  c  vector  is  too  close  to  optimum. 
The  FNR  method  does  not  reduce  the  augmented  performance 
measure  but  manages  to  salvage  the  first  iteration.   On 
iteration  number  2  the  opposite  occurs.   The  c  vector  is 
too  far  from  optimum  for  the  FNR  method  to  work.   The  MNR 
method  comes  to  the  rescue.   The  same  thing  occurs  at  the 
beginning  of  sub-problem  2.   In  these  two  sub-problems  the 
use  of  both  methods  in  combination  has  allowed  the  algorithm 
to  proceed  where  the  exclusive  use  of  either  method  by 
itself  produces  a  divergent  condition  from  which  the 
algorithm  cannot  recover. 

Figure  19  is  a  plot  of  the  performance  (final  time) 
vs.  iteration  number. 

Figures  20,  21,  and  22  are  plots  of  the  states  vs. 
time.   In  each  plot  two  curves  are  shown:   one  is  the 
expansion  for  the  state  as  computed  at  the  end  of  the  last 
sub-problem;  the  second  is  the  state  trajectory  obtained  by 
integrating  the  state  equations  with  the  optimal  control 
expansion. 
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Figure  19 

Performance  measure  vs.  iteration  number 
for  climb  performance  problem  A. 
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Figure  20 

Mach  number  vs.  time  for 
climb  performance  problem  A. 
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Altitude   vs.    time   for 
climb   performance  problem  A. 


1*  30 

01 

>Mn 

tegration        II          \ 

ai 

T3 

w  20 

—       X 

01 

^~ 

CD 

c 

■a 

5  i° 

.<o 

Q. 

/  /                 .1 

4J 

//      expansion      1  >^ 

CO 

•^   o 

1 

\        1 

//      i                      i  i 

1 

\       1 

\    2 

/ /          3                          4   » 

r— 

//      Time  (min.) 

<a 

o 

"■13  -10 

J- 

01 

Figure  22 

Vertical  flight  path  angle  vs. 
time  for  climb  performance  problem  A. 
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Close  observation  of  Figures  20,  21,  and  22  reveals  a 
trajectory  very  similar  to  that  described  in  the  introduction 
to  this  problem.   The  trajectory  begins  with  a  sub-sonic 
climb  to  an  altitude  of  33,000  feet.   The  climb  angle  during 
this  portion  of  the  climb  reaches  a  maximum  of  27  degrees. 
At  this  point  an  acceleration  is  performed  to  a  Mach  number 
of  1.53  with  the  aircraft  in  a  slight  descent.   A  "zoom" 
climb  is  then  performed  to  the  desired  altitude  of  60,000 
feet. 

Figure  23  is  a  plot  of  the  angle  of  attack  expansion 
computed  at  the  end  of  the  last  sub-problem. 
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Figure  23 

Angle  of  attack  vs.  time 
for  climb  performance  problem  A. 
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Initially,  the  angle  of  attack  is  decreasing  corresponding 
to  the  initial  acceleration  of  the  aircraft  from  a  starting 
Mach  number  of  0.6.   As  the  aircraft  rotates  to  climb 
attitude,  the  angle  of  attack  increases.   As  the  aircraft 
levels  off  for  the  supersonic  acceleration,  the  angle  of 
attack  decreases  correspondingly.   The  angle  of  attack 
begins  to  increase  again  as  the  "zoom"  climb  attitude  is 
established.   A  further  increase  is  evident  as  the  aircraft 
slows  down  in  the  climb  until  as  the  final  altitude  is 
approached,  the  aircraft  is  near  stall.   The  climb  was 
completed  in  4  minutes  and  l8vseconds. 

Figure  24  is  a  plot  of  altitude  vs.  range  where  both 
quantities  are  obtained  by  integration.   The  scales  are  not 
the  same . 
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Figure  24 

Altitude  vs.  range  for 
climb  performance  problem  A. 
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The  problem  was  also  solved  with  the  maximum  Mach 

number  restricted  to  1.0  throughout  the  flight  regime  to 

obtain  a  comparison  of  the  "zoom"  climb  technique  to  a 

totally  subsonic  climb  to  60,000  feet.   The  aircraft  was  not 

able  to  complete  the  climb.   The  altitude  of  60,000  feet  is 

apparently  above  the  service  ceiling  of  the  model.   After 

Jj  minutes  and  18  seconds,  which  was  the  time  required  to 

complete  the  climb  by  the  "zoom"  technique,  the  aircraft 

was  passing  43,000  feet  and  climbing  very  slowly. 

2.   Problem  B  -  Optimum  Climbs  to  Altitudes 
from  60,000  to  70,000  FeeT 

Table  9  depicts  the  results  for  each  sub-problem  as 
minimum  time-to-climb  profiles  are  generated  for  final 
altitudes  of  60,000  to  70,000  feet  in  thousand-foot  increments. 
For  each  sub-problem  the  results  of  the  previous  sub-problem 
were  used  as  a  first  guess  for  the  new  trajectory.   The 
stable  behavior  of  the  FNR  method  near  the  minimum  is  respon- 
sible for  the  ability  of  the  algorithm  to  generate  neighboring 
optimal  trajectories.   Thus,  in  effect,  ten  problems  were 
solved  with  minimum  effort  by  taking  advantage  of  the  results 
of  each  problem  in  turn.   If,  however,  the  change  in  end 
conditions  is  too  large,  the  MNR  method  may  be  required  to 
start  the  new  sub-problem.   This  is  the  case  in  sub-problems 
11  thru  14. 

Figure  25  is  a  plot  of  altitude  vs.  range  where 
both  quantities  are  obtained  by  numerical  integration  showing 
the  climb  trajectories  for  several  of  the  final  altitudes. 
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sub- 
problem 

number  of 

iterations 

required 

optimization 
strategy  * 

final 

altitude 

(feet) 

time  to 

climb 

(J*) 

C.P.U. 
time 

H 

3 

FPF 

60,000 

4' 18" 

2'59" 

5 

4 

FPPP 

61,000 

V2V 

3'11" 

6 

5 

FFPFF 

62,000 

k  '31'* 

3'16" 

7 

12 

FFFFFFFFFFFF 

63,000 

4 '36" 

5'02" 

8 

4 

FFFF 

64,000 

4»4l" 

3'05n 

9 

3 

MFF 

65,000 

i|»  W 

2'52" 

10 

5 

FFFFF 

66,000 

4'54" 

3'29" 

11 

4 

MFFF 

67,000 

4'  57" 

3'04" 

12 

3 

MFF 

68,000 

5'00" 

2 '50" 

13 

H 

MFFF 

69,000 

5'03" 

2 '55" 

14 

H 

MFFF 

70,000 

5' 09" 

2 '53" 

*   F  -  FNR  method 
M  -  MNR  method 


Table  9 
Results  of  climb  performance  problem  B. 
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Figure  25 

Altitude  vs.  range  for 
climb  performance  problem  B. 


The  climb  profiles  shown  in  Figure  25  reveal  the  following 
characteristics : 

a.  The  sub-sonic  climb  profiles  are  identical  for 
all  final  altitudes  for  the  initial  portion  of  the  climb; 

b.  as  the  final  altitude  increases,  the  altitude  at 
which  the  aircraft  levels  off  for  the  supersonic  acceleration 
increases  by  approximately  15  to  18  percent  of  the  final 
altitude; 
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c.  the  aircraft  performcs  a  diving  maneuver  to 
transit  the  transonic  region  with  the  maximum  dive  angle 
(13°)  reached  in  the  climb  to  60,000  feet; 

d.  as  the  final  altitude  increases,  the  aircraft 
performs  the  supersonic  acceleration  at  higher  altitudes 
with  less  altitude  loss  in  acceleration.   The  results  are 
not  in  agreement  with  standard  practice  in  which  the 
accelerations  are  generally  performed  in  level  flight  at 
the  tropopause.   The  results  are  in  agreement  with  the 
results  of  Bryson  [Ref.  2]  which  clearly  show  the  dive 
associated  with  transiting  the  transonic  region.   Bryson' s 
results  were  obtained  by  the  method  of  steepest  descent. 
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VII.   AN  AIR  COMBAT  MANEUVERING  PROBLEM 

In  this  section  the  turning  performance  of  a  supersonic 
fighter  is  considered.   First,  the  basic  aircraft  limitations 
on  maneuvering  flight  are  reviewed.   Second,  turning  perfor- 
mance in  the  horizontal  plane  is  reviewed  from  a  theoretical 
point  of  view.   The  "corner"  velocity  concept  familiar  to 

» 

fighter  pilots  is  presented.   Third,  turning  performance  in 
three-dimensions  is  discussed.   Finally,  a  three-dimensional 
problem  is  solved  in  which  a  supersonic  fighter  is  required 
to  execute  a  180°  course  reversal  in  minimum  time  with  the 
initial  and  final  altitudes  specified. 

A.   THEORETICAL  TURNING  PERFORMANCE 

1.   Aircraft  Performance  and  Maneuvering  Limitations 

A  tactical  fighter  must  be  highly  maneuverable .   An 
important  consideration  in  maneuverability  is  the  ability 
of  the  fighter  to  turn.   Two  basic  performance  criteria  in 
turning  performance  are: 

a.  radius  of  turn,  and 

b.  rate  of  turn. 

In  air  combat  situations  it  is  often  desirable  to  perform 
a  turn  so  that  the  aircraft's  radius  of  turn  (curvature)  is 
minimized  or  the  aircraft's  rate  of  turn  is  maximized.   The 
ability  of  a  fighter  to  minimize  radius  of  turn  or  maximize 
rate  of  turn  is  limited  by: 
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a.  maximum  thrust, 

b.  aerodynamic  stall,  or 

c.  maximum  allowable  load  factor. 
2.   Turns  in  the  Horizontal  Plane 

If  an  aircraft  is  restricted  to  move  in  a  horizontal 
plane  only,  turning  performance  is  easily  analyzed.   Using 
the  assumptions  given  in  Appendix  A  and  the  added  assumption 
that  v 

T  sin  a  <<  L,  (7.1) 

equations  for  lift  coefficient,  radius  of  turn,  rate  of 
turn,  and  the  thrust  required  to  maintain  level  flight  at 
constant  velocity  are  easily  derived  and  well  known.   They 
are 


2W  n 
CT  =  —^j-  ,  (7.2) 

L   pSv 


r  =  __v (7<3) 

g(n2-l)-* 


*  -*£*£.  (7.4) 


and 


T  =  — ^ v2  CD  .  (7.5) 

2  cos  a 
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where 


T  =  thrust, 

a  =  angle  of  attack, 

L  =  lift, 

C  =  lift  coefficient, 

W  =  aircraft  weight, 
e  o» 

n  =  load  factor 

p  =  air  density, 

S  =  wing  area, 

v  =  aircraft  velocity, 

g  =  gravitational  constant, 

C_  =  drag  coefficient, 

R  =  radius  of  turn,  and 

ty  =  horizontal  flight  path  angle  (heading  angle) . 


Figure  26  is  a  plot  of  lines  of  constant  thrust,  constant 
radius  of  turn,  constant  rate  of  turn,  and  maximum  lift 
coefficient  on  velocity-load  factor  CV-n)  diagrams. 


constant 
rate  of  turn  r 


/ 

/ 

"nM 

incr. 

f 

/ 

/ 

/  7x 

S 

/  / 

k/ 

/ 

//  -^ 

*  A' 

^ 

«* 

/  clm 

1 

(a) 


velocity 


constant 
radius  of  turn  R 


'    / 

7 

/     / 

1      / 

~/  Clm 

incr. 
R 

(b) 


velocity 


constant 
thrust  T 


(c) 


velocity 


Figure    26 
Velocity-Load  Factor  Diagrams 
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In  Figure  26(c)  the  constant  thrust  lines  indicate  the 

thrust  required  to  maintain  a  steady  state  turn  at  a  specific 

load  factor  and  velocity.   The  corner  velocity  v  is  defined 

as  that  velocity  at  which  an  aircraft  is  capable  of  operating 

at  maximum  lift  coefficient  CT   and  maximum  structural  load 

LM 
factor  n„  at  the  same  time.   This  is  the  velocity  which 

produces  minimum  radius  of  turn  and  maximum  rate  of  turn  as 

can  be  seen  from  Figures  26(a)  and  26(b).   The  corner 

velocity  can  only  be  maintained  in  the  steady  state  if  the 

aircraft  has  enough  thrust  available  to  allow  the  maximum 

thrust  curve  to  pass  above  the  corner  created  by  the 

CT   -  nM  boundary  intersection.   If  sufficient  thrust  is 

LM    M 
not  available  to  allow  the  corner  velocity  to  be  maintained 

in  the  steady  state,  which  is  typically  the  case,  the  air- 
craft must  either  degrade  its  turning  performance  by  moving 
off  the  boundary  intersection  until  the  maximum  thrust  curve 
is  encountered,  or  sacrifice  altitude.   In  this  case  the 
velocity  for  maximum  rate  of  turn  is  larger  than  the  velocity 
for  minimum  radius  of.  turn. 

As  can  be  seen  from  the  previous  discussion, 
optimization  techniques  are  not  required  to  analyze  turns  in  the 
horizontal  plane.   This  arena,  however,  is  an  excellent 
place  to  test  an  optimization  technique  which  is  being 
considered  for  use  in  solving  more  complicated  problems 
involving  three-dimensional  maneuvering.   With  this  in  mind, 
two  problems  involving  turns  in  the  horizontal  plane  were 
solved  by  the  epsilon  method  and  the  answers  compared  to 
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theoretical  results.   In  one  problem  an  aircraft  was 
required  to  perform  a  horizontal  turn  with  minimum  radius. 
In  a  second  problem  the  aircraft  was  required  to  turn 
through  a  given  heading  change  in  minimum  time  which  is 
equivalent  to  maximizing  rate  of  turn.   The  aircraft  was 
given  a  large  thrust  capability  so  that  the  turns  were  not 
thrust  limited.   The  mathematical  model  used  in  the  problems 
is  given  in  Appendix  A.   The  model  is  an  accurate  three- 
degree-of-freedom  model  of  the  aircraft's  motion  with 
maneuvering  limitations  included.   Prom  the  previous 
discussion  the  aircraft  should  have  performed  the  turn 
in  both  cases  at  the  corner  velocity  where 


c  *   pS(Cn  tana  +  CT  ) 
L     DM    S    LM 


(7.6) 


In  this  equation  a  is  the  angle  of  attack  for  maximum  lift 
coefficient  and  Cn  is  the  drag  coefficient  for  maximum  lift 
coefficient.   The  epsilon  method  solved  both  problems 
successfully.   In  each  case  the  optimum  trajectory  involved 
a  turn  at  the  corner  velocity.   Thus,  the  ability  of  the 
second-order  epsilon  method  to  handle  problems  of  this  type 
was  demonstrated. 

3.   Turns  in  Three  Dimensions 

The  analysis  of  turning  performance  in  three 
dimensions  is  quite  complicated.   In  this  regime  modern 
optimization  techniques  are  the  only  method  of  solving 


118 


meaningful  problems.   Even  optimization  techniques  are  apt 
to  have  a  difficult  time  with  three-dimensional  maneuvering 
problems  using  realistic  models  of  the  aircraft's  motion 
because  of  large  computer  time  and  storage  requirements.   In 
this  section  the  epsilon  method  is  used  to  solve  an  impor- 
tant three-dimensional  maneuvering  problem  often  encountered 
in  air  combat. 

In  many  air  combat  situations  a  fighter  pilot  is 
faced  with  the  requirement  to  reverse  his  course  as  rapidly 
as  possible.   Generally,  the  pilot  has  in  mind  a  specific 
altitude  at  which  he  would  like  to  complete  the  maneuver 
which  is  dictated  by  his  desire  to  track  an  enemy  aircraft 
or  perform  some  attacking  maneuver.   With  this  in  mind,  the 
problem  posed  for  solution  by  the  epsilon  method  involves  a 
supersonic  fighter  which  is  required  to  execute  a  180° 
reversal  in  minimum  time.   The  aircraft  must  begin  the 
maneuver  in  level  flight  and  recover  in  level  flight  at  the 
entry  altitude.   The  accepted  maneuvers  used  to  accomplish 
this  task  developed  over  years  of  combat  experience  are  the 
high-speed  yo-yo  and  the  low-speed  yo-yo  maneuvers.   If  the 
aircraft  begins  a  reversal  at  a  flight  speed  higher  than 
its  corner  velocity  a  high-speed  yo-yo  is  called  for  and 
vice-versa.   A  high  speed  yo-yo  consists  of  a  climbing  turn 
followed  by  a  descending  turn.   A  low  speed  yo-yo  consists 
of  a  descending  turn  followed  by  a  climbing  turn.   If  the 
aircraft  begins  a  reversal  at  its  corner  velocity,  a  level 
turn  is  called  for.   The  purpose  of  applying  the  epsilon 
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method  to  this  problem  is  to  either  confirm  or  challenge 
the  effectiveness  of  these  experimentally  developed 
maneuvers  by  the  use  of  an  optimization  technique.   The 
assumptions  applied  to  the  problem,  the  coordinate  system 
and  nomenclature  used,  and  the  derivation  of  the  equations 
of  motion  are  presented  in  Appendix  A. 

B.   PROBLEM  FORMULATION 
1.   State  Equations 

The  state  equations  derived  in  Appendix  A  are 


m  =  wir  Thcosct  -  f  slnT  -  ffr^  m2°d  (7.7) 

e  e 

^gnsincf,  (   8) 

y   a  Mcosy 


'    £  ncos(j)  _  g  cosy  ,y  q) 

T   a   M     a   M  \i*yj 


h   =  fp.   siny   .  (7.10) 

HL 


The  states  are  Mach  number  M,  horizontal  flight  path  angle 
i|>,  vertical  flight  path  angle  y,  and  normalized  altitude  h. 
The  controls  are  bank  angle  4>,  normalized  thrust  Th,  angle 
of  attack  a,  and  load  factor  n. 

In  addition  to  the  four  state  equations,  an  equality 
constraint  which  must  be  satisfied  is 


0  =  ~£L  Thsina  +  ^f-  M^CL  -  n  (7.11) 

e  e 
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This  equation  is  a  result  of  the  definition  of  load  factor 
n  and  is  derived  in  Appendix  A.   It  is  possible  to  use 
equation  (7.11)  to  eliminate  the  load  factor  control  from 
the  state  equations,  but  this  is  not  desirable  for  two 
reasons:   first,  the  resulting  state  equations  would  be 
further  complicated  thus  increasing  the  analytic  workload 
required  to  compute  first  and  second  partial  derivatives; 
second,  the  incorporation  of  the  important  load  factor 
inequality  constraint  would  be  unnecessarily  complicated. 

Parameters  considered  constant  for  this  problem  are 
the  gravitational  constant  g,  aircraft  maximum  thrust  ThM, 
aircraft  weight  W  ,  speed  of  sound  a,  base  altitude  Hj- ,  air 
density  p,  and  aircraft  wing  area  S.   These  constants  are 

g   =  32,1725  ft. /sec.2 

ThM  =  21,000   lbs. 

W   =  39,000  lbs. 
e 

a   -  1077.8  ft. /sec.  (7.12) 


HL  =  10,000  ft. 
p   =  0.001756  slugs/ft.3 
S   =  400  ft.2 
It  is  convenient  to  define  the  constants 

c-  =  f  (7.13) 

1   a 


A  gThM 

c2  "  W  a 
e 


(7.14) 


c  =  gPSa  .  (7.15) 

c3   2W  w.-io, 

e 


121 


1   \ 


oil 

A 

We 

C5 

A 

pSa2 

2W 
e 

(7.16) 


(7.17) 


With  these  simplifying  definitions  incorporated  the  state 
equations  are 


*  ? 

M  =  c2Thcosa  -  c-^iny  -  c-M  CD  (7.18) 


^=cl^  (7.19) 


ncos(j>  _    cosy  ,   2  . 

i   cl   M     Cl   M  w.^u; 


h  =  |H  siny  (7.21) 

HL 


and  the   additional  equality   constraint   is 


0   =   CjjThsina  +   c5M2CL  -  n      .  (7.22) 


2.   Tabular  Functions 

Tables  are  used  for  lift  and  drag  coefficient  as 
functions  of  Mach  number  and  angle  of  attack  for  a  typical 
supersonic  fighter.   Excerpts  from  these  tables  are  provided 
in  Appendix  B. 
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3.   Performance  Measure 

The  performance  measure  for  this  problem  is 

T 


- f    dt.  (7.23) 


J 


k.      Inequality  Constraints 

The  controls  must  satisfy 


0  <  Th  <  1  (7.24) 


0  i  a  <  aM  =  24°  (7.25) 


°lninM=6'5g's  (7.26) 

0  <  <J>  <  u  .  (7.27) 

The  minimum  allowable  normalized  thrust,  angle  of  attack, 
and  load  factor  are  approximated  by  zero  as  these  constraints 
are  not  anticipated  to  be  active.   A  zero  value  of  the 
lower  bound  simplifies  the  associated  penalty  term  in  the 
augmented  performance  measure.   The  maximum  angle  of  attack 
is  fixed  at  a  value  slightly  higher  than  the  angle  of  attack 
for  maximum  lift  coefficient  as  given  in  the  tabular  data 
thus  simulating  aerodynamic  stall.   The  structural  load 
factor  upper  bound  is  6.5  g's,  a  standard  value  from  fighter 
aircraft  operational  limitations.   The  bank  angle  constraints 
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were  required  to  keep  the  algorithm  from  generating  bank 
angles  greater  than  180°. 
5.   End  Conditions 

The  initial  conditions  are 

M(0)  =  0.9  (7.28) 

*(0)  =  0  (7.29) 

Y(0)  =  0  (7.30) 


h(0)  =  hQ  =  15,000  ft.  (7.3D 


The  final  conditions  are 


*-(T)  =  tt  (7.32) 


Y(T)  =  0  (7.33) 


h(T)  =  hQ   .  (7.34) 


C.   THE  EPSILON  METHOD  FORMULATION 

1 .   The  Augmented  Performance  Measure 

Using  the  penalty  functionals  described  in  Section 
III  for  the  inequality  constraints,  the  augmented  perform- 
ance measure  is 
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.T 
J_   =  I       dt 


v0 


T 
+  -y       ("m  -   c2Thcosa  +   c   siny  +   c-M2CD]2   dt 


+  |y       Tc^Thsina  +   c5M2CL  -  nl  2   dt 


(7.35) 


•T 

Si 


where  e  and  r  are  weighting  factors,  and  6  is  a  constant 
used  to  make  minor  adjustments  to  the  admissible  regions. 
The  required  elements  of  w  are 

wk  =  [\  ~  ('2I\C0S\  +  cisin\  +  c^k\  ]  [itT  »  k=1'2»  • ' '  »K  (  7 '  36 ) 


=  [^-cl^s%][f]%   '  k=1'2 K  (7<37) 
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w 


w 


2K+k 


3K+k 


=    [\  '  Cl 


a  cos  4  cosy 


\ 


\ 


-][~]%.  k=l,2,...,K  (7.38) 


■  [\  -  r sin\]  [t f   >k=1^-->K 


w 


4K+k 


w 


-r^-i 


IK 


5K+k        |_"~S 


'(rAt)'*      ,      k=l,2,...,K 


w 


6K+k 


_  r2ak 


IK 


w 


7K+k 


w 


8K+k 


M 

2*k      - 
IT  ~  1 


K 


'(rAt)2      ,      k=l,2,...,K 


(rAt)**      ,      k=l,2,...,K 


P(rAt)2      ,    k=l,2,...,K 


(7.39) 


=    [ci,T\slnak+  c5\\    ~\]  [tP    ,k=l,2,...,K  (7.40) 


(7.4l) 
(7.42) 

(7.^3) 
(7.44) 


% 


W9K+1   =    [(K-l)Atr 


(7.45) 


2.   Functional  Expansions,  Unknowns,  and 
Partial  Derivatives 

The  state  and  control  expansions  are  of  the  same 

form  as  in  previous  problems  and  are  not  shown.   The  elements 

of  the  c  vector  are 


T 
c     =    (a,  ta.0 , .  .  .  ,a    , b,  , o0 , . . .  ,b    , c,  , c0 ,  . .  . , cM, d,  , d5 ,  . .  .  , dM, 


ll»~2 


lm»    1»    2 


m»"l»v'2- 


'M»"l»"2 


W 


ei »eo  > •• • »eM»i »*o» • • •»  ^m»St i So  » • • • » Sm>"t  »"o » • • • »"m» 


1,"21 


*M»    1»    2 


j^»&1»&2 


=M»    1»    2- 


M 


K»4'1»<<'K,Th1,ThK>a1,aK,n1,nK,At) 


W- 


(7.46) 
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where  the  a  '  s  represent  the  Mach  number  expansion  coeffi- 
cients, the  b  's  represent  the  horizontal  flight  path 
coefficients,  the  c  '  s  represent  the  vertical  flight  path 
angle  coefficients,  the  d  '  s  represent  the  altitude  coeffi- 
cients, the  e  '  s  represent  the  bank  angle  coefficients,  the 
f  's  represent  the  thrust  coefficients,  the  8L.'s  represent 

the  angle  of  attack  coefficients,  and  the  h  's  represent 

'  m 

the  load  factor  coefficients. 

The  first  and  second  partial  derivatives  of  the 
tabular  functions  for  lift  coefficient  and  drag  coefficient 
with  respect  to  their  independent  variables  and  equations 
(7.36)  thru  (7.^5)  with  respect  to  the  c  vector  are  obtained 
in  the  same  manner  as  in  previous  problems. 


D.   RESULTS 

Three  problems  were  solved.   In  problem  A  the  aircraft 
must  perform  the  180°  reversal  in  minimum  time  starting 
from  an  initial  Mach  number  of  0.9.   In  problem  B  the  air- 
craft must  perform  the  reversal  starting  from  its  corner 
Mach  number  which  from  equation  (7.6)  is  0.708.   In  problem 
C  the  aircraft  must  perform  the  reversal  starting  from  an 
initial  Mach  number  of  0.5.   In  all  problems  8  coefficients 
for  each  expansion  and  21  time  points  (20  time  intervals) 
are  used. 

1.   Problem  A 

Since  the  initial  Mach  number  is  above  the  corner 
Mach  number,  a  high-speed  yo-yo  maneuver  is  called  for  by 
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accepted  tactics.   With  this  in  mind  an  initial  guess  of 
the  c  vector  was  made  to  reflect  this  type  of  maneuver. 
Accordingly,  the  following  coefficients  were  given  non-zero 
initial  values: 

a1   =  -0.309 

c,  -  0.524 

d,  =  0.333 

1  (7.^7) 

e1  =  1.047 

g1  =   0.262 

h±   =   3.000 

The  remaining  initial  values  for  the  c  vector  were: 

remaining  expansion  coefficients  =  0 

M(T)  =  0.9 

4)(0)  =  0° 

<f>(T)  =  0° 

Th(0)  =  0.88  (30,000  lbs.) 

Th(T)  =  0.88  (30,000  lbs.)  (7.48) 

a(0)  =  5° 

a(T)  =  5° 

n(0)  =  1.0 

n(T)  =  1.0 

At  =  12/20  sec.  (T  =  12  sec.) 

The  problem  was  solved  in  six  sub-problems.   It  took 
17.65  seconds  to  complete  the  turn.   Figures  27  thru  30 
are  plots  of  the  control  expansions  as  computed  at  the  end 
of  the  last  sub-problem. 

From  Figure  28  it  is  seen  that  the  maximum  thrust 
constraint  is  active  for  the  first  10  seconds  of  the  turn. 
At  t  k  6  seconds  there  is  a  short  period  in  which  the 
thrust  is  slightly  inadmissible.   This  is  due  to  the  use 
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Figure  27 

Bank  angle  vs.  time  for 
turning  performance  problem  A. 
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Figure  28 

Thrust  vs.  time  for 

turning  performance  problem  A. 
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Figure  29 

Angle  of  attack  vs.  time 

for  turning  performance  problem  A. 
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Figure  30 

Load  factor  vs.  time  for 
turning  performance  problem  A. 
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of  the  factor  8   =  1.03  in  the  inequality  constraint  penalty 
terms  in  equation  (7.35)  which  has  the  effect  of  slightly 
increasing  the  size  of  the  admissible  region.   This  is 
desirable,  however,  as  the  epsilon  method  generates  only  an 
approximation  to  the  optimal  control  from  which  the  true 
optimal  control  must  be  deduced.   It  is  easier  to  recognize 
an  optimal  control  expansion  which  is  on  a  constraint 
boundary  with  the  factor  <S  included.   As  shown  in  Figure  6 
in  Section  V,  6  =  1,03  is  the  proper  choice  for  a  final 
K  =30.   During  the  last  portion  of  the  turn,  the  aircraft 
is  operated  at  the  angle  of  attack  for  maximum  lift  coeffi- 
cient (approximately  22°).   The  bank  angle  and  load  factor 
constraints  are  not  active  during  the  maneuver. 

Figures  31  thru  3^  are  plots  of  the  states  vs.  time. 
In  each  plot  two  curves  are  shown:   one  is  the  expansion  for 
the  states  as  computed  at  the  end  of  the  last  sub-problem; 
the  second  is  the  state  trajectory  obtained  by  numerically 
integrating  the  state  equations  with  the  optimal  control 
expansions.   An  observation  of  these  plots  reveals  that  a 
high-speed  yo-yo  maneuver  is  performed  although  the  altitude 
excursions  are  not  as  large  as  this  author  expected.   The 
optimization  procedure  settles  on  a  nearly  level  hard  turn 
at  high  load  factors,  steep  bank  angles,  and  maximum  thrust 
for  the  majority  of  the  turn.   When  the  maximum  thrust 
boundary  is  not  active,  the  aircraft  flies  at  the  angle  of 
attack  for  maximum  lift  coefficient. 
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Figure    31 

Mach  number  vs.    time 

for  turning  performance  problem  A. 
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Figure    32 

Horizontal   flight  path  angle   vs.    time 
for  turning  performance   problem  A. 
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Figure    33 

Vertical  flight  path  angle  vs.  time 
for  turning  performance  problem  A. 
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Figure  3^ 

Altitude  vs.  time  for 

turning  performance  problem  A. 


2.   Problems  B  and  C 

Figures  35,  36,  and  37  are  plots  of  cross-range  vs. 
range,  altitude  vs.  cross-range,  and  altitude  vs.  range 
obtained  by  integrating  the  equations  of  motion  with  the 
optimal  control  expansions  found  in  problems  A,  B,  and  C. 
An  observation  of  Figures  35,  36,  and  37  reveals  that  the 
expected  maneuvers  are  performed  for  each  initial  Mach 
number.   In  problem  B  the  aircraft  performs  an  essentially 
level  turn  from  an  initial  Mach  number  equal  to  its  corner 
Mach  number  at  this  altitude.   In  problem  C  the  aircraft 
performs  a  low-speed  yo-yo  maneuver  losing  a  maximum  of 
800  feet  after  90°  of  turn  from  an  initial  Mach  number 
below  its  corner  Mach  number.   In  problems  A  and  C,  however, 
the  aircraft  does  not  go  through  as  much  of  an  altitude 
excursion  as  anticipated  by  the  author.   Since  in  fighter 
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Range  (ft.  x  10 r) 

2  4  6 


8 


problem  B 


problem  A 


Figure   35 

Cross-range   vs.    range   for 

turning  performance  problems   A,    B,    and  C. 
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Altitude  vs.  cross-range  for 

turning  performance  problems  A,  B,  and  C. 
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Altitude  vs.  range  for 

turning  performance  problem  A,  B,  and  C- 
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tactics,  however,  there  are  no  rules  on  the  amount  of 
altitude  which  should  be  gained  or  lost  in  a  yo-yo  maneuver, 
a  quantitative  evaluation  of  the  results  is  purely  subjec- 
tive.  The  important  result  is  that  the  optimization  tech- 
nique did  require  the  aircraft  to  perform  the  high-speed 
and  low-speed  yo-yo  maneuvers  predicted  by  accepted  tactics. 
The  accepted  tactics  are,  therefore,  qualitatively  correct. 
The  optimal  times  required  to  complete  the  turn  for 
each  of  the  three  problems  are  given  in  Table  10. 


Optimal  time 
for  reversal 

problem  A 
problem  B 
problem  C 

17. 6  sec . 

20.8  sec. 

24.9  sec. 

Table  10 
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VIII.   SUMMARY  AND  CONCLUSIONS 

A  number  of  realistic  problems  in  aircraft  and  missile 
performance  optimization  have  been  solved  by  the  use  of  a 
second-order  epsilon  method.   The  mathematical  models  have 
portrayed  aircraft  and  missile  dynamics  In  an  accurate 
manner  with  particular  emphasis  placed  on  the  modeling  of 
performance  and  maneuvering  limitations. 

The  state  and  control  inequality  constraints  generated 
by  these  limitations  have  been  handled  by  a  new  computa- 
tionally superior  penalty  functional.   Three  desirable 
theoretical  properties  of  this  penalty  functional  have  been 
shown. 

A  full  Newton-Raphson  method  for  minimizing  the  aug- 
mented performance  measure  has  been  shown  to  be  computation- 
ally feasible  and  superior  in  certain  situations  to  the 
"modified"  Newton-Raphson  method  proposed  elsewhere. 

The  following  observations  are  significant  with  respect 
to  the  second-order  epsilon  method. 

a.  The  MNR  method  is  relatively  insensitive  to  the 
starting  values  of  the  unknowns  c.  The  FNR  method  diverges 
for  starting  values  of  c  which  are  far  from  optimum. 

b.  Once  c  is  close  to  optimum,  the  FNR  method 
converges  rapidly  whereas  typically  the  MNR  method  either 
diverges,  oscillates,  or  converges  slowly  at  best. 

c.  In  relatively  simple  problems  the  MNR  method  is 
capable  of  obtaining  a  solution  by  itself.   In  more  difficult 
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problems  such  as  those  solved  in  this  dissertation,  a 
combination  of  the  two  methods  is  required.   Typically,  the 
most  effective  procedure  involves  using  the  MNR  method 
initially  followed  by  the  FNR  method  when  successive  itera- 
tions yield  "small"  improvements  in  the  augmented  perform- 
ance measure.   In  other  rare  cases  where  the  initial  guess 
for  the  c  vector  is  close  to  optimum,  the  FNR  method  must 
be  used  initially. 

d.   The  power  of  the  FNR  method  close  to  the  minimum 
can  be  used  to  advantage  to  obtain  with  minimum  effort 
optimal  control  and  state  trajectories  for  problems  with 
neighboring  end  conditions  by  using  the  solution  to  a  basic 
problem  as  a  first  guess  for  the  new  problem. 

The  solutions  to  the  problems  solved  have  a  number  of 
applications.   In  the  missile  intercept  problem  (Section  V) 
minimum-time  optimal  trajectories  may  be  used  as  a  basis  for 
comparison  with  the  performance  of  more  practical  sub- 
optimal  controllers  such  as  proportional  navigation  for  a 
short  range  air-to-air  missile .   In  the  three-dimensional 
turn-reversal  problem  the  qualitative  optimality  of  an 
experimentally  developed  air  combat  maneuver  is  shown  for 
the  first  time.   A  significant  lesson  to  be  learned  from 
the  results  of  this  problem  is  the' importance  of  thrust  in 
comparison  to  lift  coefficient  in  the  maneuvering  capability 
of  a  fighter  aircraft.   Thus,  an  optimization  method  of  the 
type  used  in  this  work  applied  to  realistic  performance 
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problems  has  application  in  the  evaluation  of  tradeoffs  in 
the  design  of  future  flight  vehicles. 
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APPENDIX  A 
MATHEMATICAL  MODELS 

In  this  Appendix  the  mathematical  models  used  in  the 

problems  are  derived.   In  Section  A.l  the  basic  equations 

of  motion  of  an  aircraft  in  three  dimensions  are  derived 

under  appropriate  assumptions.   This  model  is  used  in  the 

problem  solved  in  Section  VII.   In  Section  A. 2  the  aircraft 

is  restricted  to  move  in  the  horizontal  plane  only  and  the 

appropriate  adjustments  are  made  to  the  three-dimensional 

model.   In  Section  A. 3  the  aircraft  is  restricted  to  move 

in  the  vertical  plane  only  and  the  appropriate  adjustments 

are  made  to  the  three-dimensional  model.   This  model  is 

used  in  the  problem  solved  in  Section  VI.   In  Section  A. 4 

the  mathematical  model  for  the  missile  intercept  problem 

solved  in  Section  V  is  derived. 

1.   The  Mathematical  Model  for  an  Aircraft  Maneuvering 
in  Three  Dimensions 

In  this  section  the  basic  three-degree-of-freedom 

equations  of  motion  of  an  aircraft  are  derived.   The 

assumptions  are 

a.  the  earth  is  flat, 

b.  the  aircraft  is  a  point  mass, 

c.  the  mass  of  the  aircraft  is  a  constant, 

d.  the  aircraft  is  always  in  balanced  flight, 

e.  the  aircraft  rolls  about  its  velocity  vector, 
and     f.   acceleration  due  to  gravity  is  a  constant. 


Ill 


»   \ 


The  coordinate  system  and  notation  are  presented  below, 


Figure  39 
Aircraft  Coordinate  System 


Three  axis  systems. are  drawn  in  Figure  39.   They  are 


a. 

b. 


(X,  Y,  H)     a  fixed  inertial  axis  system; 


(x- 


Y' 


c.   (x,  y,  z) 


Z1)   a  non-rotating  axis  system  fixed 
to  the  center  of  mass  of  the 
aircraft; 

a  rotating  axis  system  fixed  to 
the  center  of  mass  of  the  aircraft; 
the  x  axis  is  oriented  in  the 
direction  of  the  aircraft's  velocity 
vector;  the  y  axis  points  out  the 
right  wing. 
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The  attitude  of  the  aircraft  is  given  by  four  angles  as 
follows: 

a.  a  rotation  \\>   in  the  X'Y'  (horizontal)  plane; 

b.  a  rotation  y   in  the  xZ'  (vertical)  plane; 

c.  a  rotation  <J>  in  the  yz  plane; 

d.  a  rotation  a  in  the  xz  plane. 

The  three  angles  ty,    ys    and  4>  are  the  Euler  angles  (Ref.  19). 
The  angle  a  is  the  angle  of  attack  of  the  aircraft  using 
the  thrust  line  as  a  reference.   The  remaining  notation  is 

a.  forces; 

L  =  lift, 

D  =  drag, 

T  =  thrust, 

Th  =  normalized  thrust, 

W  =  gross  weight, 

b.  angles; 

a   =  angle  of  attack  of  the  thrust  line 
measured  in  the  xz  plane, 

ct_   =  angle  of  attack  for  maximum  lift  coefficient, 

Y   =  vertical  flight  path  angle  measured  in  the 
xZ'  plane, 

(j>   =  bank  angle  measured  in  the  yz  plane, 

ty       =  horizontal  flight  path  angle  measured  in 
the  X'Y'  plane, 

c.  rates; 

p   =  roll  rate  measured  in  the  yz  plane, 
q   =  pitch  rate  measured  in  the  xz  plane, 
r   =  yaw  rate  measured  in  the  xy  plane, 
w   =  angular  velocity  of  the  xyz  system  with 
respect  to  the  X'Y'Z1  system, 

d.  other  parameters; 

v   =  velocity, 

m   =  mass, 

g   =  gravitational  constant, 


1^3 


p   =  air  density, 

S   =  aircraft  wing  area, 

C   =  lift  coefficient, 

CD  =  drag  coefficient, 

H   =  altitude, 

h   =  normalized  altitude, 

R   =  radius  of  turn, 

M   =  Mach  number, 

M   =  Corner  Mach  number, 
c  ' 

a   =  speed  of  sound, 

v   =  corner  velocity, 
c 

n   =  load  factor, 

e   =  unit  vector  (with  appropriate  subscript 

indicating  direction), 

e.   subscripts; 

M   =  maximum  value , 
L   =  minimum  value . 

The  equations  of  motion  are  derived  following  the 

methods  used  in  Reference  19-   Starting  with  Newton's 

second  law 


dv 

?  =  m  df  =  m 


5v 
6t 


+  to  X  V 


(A.l) 


where  -r-r-  is  defined  as  the  time  derivative  in  the  xyz  system, 
The  aircraft  velocity  and  acceleration  may  be  written 


v  =  ve 


(A. 2) 


dt  =  I   =   v?x  +  v?x  > 


(A. 3) 


where 


1^4 


v?x  = 


6v 

IF 


(a.  4) 


and 


ve  =  w  x  v  . 
~x   ~ 


(A. 5) 


The  angular  velocity  of  the  xyz  system  with  respect  to  the 
non-rotating  frame  X'Y'Z1  is  given  by 


w  =  Pe  +  qe   +  re   . 


(A. 6) 


At  this  point,  relations  between  the  angular  rates  p,  q, 
and  r  and  the  angular  rates  of  change  of  the  Euler  angles 
are  required.   These  relations  are  purely  trigonometric  in 
nature  and  are  derived  in  Reference  19.   In  matrix  notation 
they  are 


P 

q 

r 


1  0  -siny 

0        cos<t>      sin<{)    cosy 

0       -sin<{>     cos<{)    cosy 


Y 


(A. 7) 


Substituting  equations  (A. 7)  into  equation  (A. 6),  we  obtain 


no  =    (<J>   -   ij)siny)e      +    (ycos<t>   +   ij>sin<}>   cosy)e 
~  ~x  ~y 


+    (ij>cos<{>   cosy   -  ysin<{>)e    . 


(A. 8) 
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Using  equations  (A. 8)  and  (A. 2)  the  product 

w  x  v  =    (ij>cos<J>   cosy  -  ysin<J>)ve      -    (ijjsin<{>   cosy  +  Ycos4>)ve 
~       —  ~y  - z 

(A. 9) 
is  formed.   Isolating  thrust  and  weight  components  in  the 
xyz  system, 


T    =   Tcosa  ,  (A. 10) 


Ty    =   0  ,  (A. 11) 

T    =   -Tsina  ,  (A. 12) 


W"    =   -Wsiny  ,  (A. 13) 

ex       e 


W  =   W  cosy  sin<J)  ,  (A.14) 

ey 

W  =  W  cosy  cosiji  ,  (A. 15) 

ez      e 


equation  (A.l)  may  be  written  in  component  form  as 


Tcosa  -  D  -  W  siny       =  mv  ,  (A.16) 

W  cosy  sintf)  =  mv(\|/cos<{)  cosy  -  Ysin<J>) ,  (A. 17) 

W  cosy  cost))  -  Tsina  -  L  =  -mv(^sinij)  cosy  +  ycos<i>).  (A.  18) 


The  equations  (A.l6)  thru  (A.l8)  are  the  basic  equations  of 
motion.   To  apply  optimization  methods,  It  is  desirable  to 
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transform  these  equations  into  state  variable  format. 
First,  lift  and  drag  coefficients  are  defined  by  the 
expressions 


L  =  CTJgpv2S  ,  (A. 19) 


D  =  C^pv^S  .  (A. 20) 


Second,  it  is  convenient  to  introduce  the  normalizing 
expressions 

v  =  aM  ,  (A.  21) 


T  =  ThMTh  .  (A. 22) 


Substituting  equations  (A. 19)  thru  (A. 22)  into  equations 
(A.16)  thru  (A.18)  along  with  the  expression 


We  =  mg  ,  (A. 23) 


we  obtain 


3  W 

ThMThcosa  -%CDp(aM)*S  -  w^siny  =  -H.  aM  (A.  24) 

W 

W  cosysintf)   =  —  aM(ipeos(J>   cosy  -  ysin<f>)  (A. 25) 

e  g    - 

W   cosy   cos<f>   -   ThMThsina   -  %C-.p(aM)    S  (A. 26) 

W 
e        /  •  \ 

=  —  aM(iJjsin$    cosy   +   ycos<{>). 
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Solving  equations  (A. 24)  thru  (A. 26)  for  M,  i> ,  and  y   yields 
the  state  variable  format 


M  =  -  Thcosa  -  £  siny  -  ^^   M2Cn  (A. 27) 

W  a  a        2We    D 

e 

1  _  gThM  Thsina  sin*   ,  gpSa  MCLsln(t>  ,.  p|n 

v    Wea    Mcosy         2We     cosy  K*.*o) 

gTh 

•  =    M  Thsina  cos^    gpSa         _  g_  cosy          ,     . 

Y    Wa     M  2V/  U0LC0S(P   a  M            ^A.^y; 

e  e 


In  addition  the  position  of  the  aircraft  (center  of  mass 
location)  may  be  required  from  some  fixed  reference  point. 
To  this  end  three  additional  state  equations  are 

X  =  aMcosy  cos^  ,  (A. 30) 

Y  =  -aMcosy  sinij/  ,  (A. 3D 

H  =  aMsiny  .  (A. 32) 

It  is  convenient,  also,  to  define  the  load  factor  n  as 


-  TOTAL  LIFTING  FORCE  ,.    „, 

n        WEIGHT  W'33) 


L  +  Tsina 
We 


(A. 3^) 


ThM  ~  2   0 

_M  Thsina  +  §^-M2CL  (A.  35) 
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With  equations  (A. 35)  incorporated  into  equations  (A. 27) 
thru  (A. 29)  the  state  equations  are 


M  =  -^-f  Thcosa  -  fsiny  -  ^^  M2CD          (A. 36) 

e  e 

•   £  nsln$  ,    , 

w        a   Mcosy  KA.31J 


•   £  ncosft    £  cosy  ,    8) 

'    a   M      a   M  ^a.jo; 


The  mathematical  model  for  the  three-dimensional 
reversal  problem  solved  in  Section  VII  includes  the  state 
equations  (A. 36),  (A. 37),  (A. 38),  and  (A. 32).   In  addition, 
equation  (A. 35)  must  be  satisfied.   This  equation  is  written 
in  the  form 


ThM  oSa2   2 

0  =  -^4  Thsina  +  ^_  M*c  _  n  .  (A.  39) 

e  e 


The  states  are  Mach  number  M,  horizontal  flight  path  angle 

i|>,  and  vertical  flight  path  angle  y.   The  controls  are 
bank  angle  <j>,  normalized  thrust  Th,  angle  of  attack  a, 
and  load  factor  n.   The  purpose  of  introducing  load  factor 
as  an  independent  control  through  equation  (A. 35)  vice  using 
the  state  equations  (A. 27)  thru  (A. 29)  is  to  simplify  the 
state  equations  and  the  incorporation  of  the  structural  load 
factor  constraint. 
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The  following  Inequality  constraints  are  imposed 
on  the  controls: 

0  <  Th  <  1  ,  (A.40) 


0  <  a  <  ar  (A. Hi) 


0  i  n  1  nM  •  (A. 42) 


The  lift  and  drag  coefficients  are  given  as  tabular 

functions  of  Mach  number  and  angle  of  attack.   Reynold's 

number  effects  are  neglected.   The  parameters  considered 

constant  for  the  problem  in  Section  VII  are  the  gravitational 

constant  g,  maximum  thrust  ThM,  aircraft  gross  weight  W  , 

the  speed  of  sound  a,  air  density  p,  and  wind  area  S. 

2.   The  Mathematical  Model  for  an  Aircraft  Maneuvering 
in  the  Horizontal  Plane 

In  this  section  the  state  equations  (A. 36)  thru 

(A. 38)  are  applied  to  an  aircraft  restricted  to  maneuver 

in  a  horizontal  plane.   The  appropriate  assumptions  are 

Y  =  0  (A. 43) 

Y  =  0  (A. 44) 
H  =  H  (A. 45) 
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Applying  equations  (A. 43)  thru  (A. 45)  to  equation  (A. 38), 
we  obtain 


n  ■  TT^T  •  (A. 46) 

COS(j> 


Substituting  equations  (A. 43)  and  (A. 46)  into  equations 
(A. 36)  and  (A. 37),  we  may  write  the  state  equations  as 


M  =  ^TiT  Thcosa  "  BzTL  m2°D   »  U'h7) 

e  e 


i  =  ^     •  (A. 48) 


The  mathematical  model  for  the  two  dimensional 
minimum  time  and  minimum  radius  of  turn  problems  referred 
to  in  Section  VII  includes  the  state  equations  (A. 47)  and 
(A. 48).   In  addition,  equation  (A. 35)  must  be  satisfied. 
Using  equation  (A. 46),  equation  (A. 35)  is  written  in  the 
form 


Th  2 

0  =  -—•  Thsina  cos*  +  ^-  CLM2cos4>  -  1  .     (A. 49) 


The  states  are  Mach  number  M,  and  horizontal  flight  path 
angle  i>.      The  controls  are  bank  angle  4>,  angle  of  attack  a, 
and  thrust  Th.   It  is  possible  to  eliminate  one  control  by 
substituting  equation  (A. 49)  into  equation  (A. 48).   The  use 
of  equation  (A. 49)  as  an  additional  equality  constraint  is 
preferred,  however,  because  the  state  equations  are  simpler 
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and  the  required  control  inequality  constraints  are  simpler 
to  incorporate. 

The  following  inequality  constraints  are  imposed  on 
the  controls: 

0  <  Th  <  1   ,  (A. 50) 

0  <  *  <  <t>M  =   cos'1[^]  ,  (A. 51) 


0  1  a  1  aM   '  (A. 52) 


The  lift  and  drag  coefficients  are  given  as  tabular 

functions  of  Mach  number  and  angle  of  attack.   The  parameters 

considered  constant  for  the  problem  are  the  same  as  those 

listed  for  the  three  dimensional  model  described  in  Section 

A.l. 

3.   The  Mathematical  Model  for  an  Aircraft 
Maneuvering  in  the  Vertical  Plane 

In  this  section  the  state  equations  (A. 27)  thru  (A. 29) 

are  applied  to  an  aircraft  restricted  to  maneuver  in  the 

vertical  plane.   The  appropriate  assumptions  are 

*  =  0  (A. 53) 

and 

i  -  0.  (A. 54) 
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Substituting  equations  (A. 53)  and  (A. 5*0  into  equations 
(A. 28)  and  (A. 29),  we  may  write  the  state  equations  as 

M  =  fl—E  Thcosa  -  §  siny  -  fg^  M2CD  (A. 55) 

e  e 

KTh 
•  =  ,-  M  Thsina   gpSa  MC   _  £  cosy  ,    6) 

Y   W  a     M     2W   '  °L   a   M  lA°o; 

e  e 

The  mathematical  model  for  the  two  dimensional  climb 
performance  problem  solved  in  Section  VI  includes  the  state 
equations  (A. 55)  and  (A. 56)  along  with  state  equation  (A. 32) 
It  is  convenient,  however,  to  introduce  the  following 
relations  into  the  state  equations: 


(A. 57) 


a 

=  _£_ 

po 

Th 

ThM 
We 

h 

H 
"HL 

(A. 58) 


(A. 59) 


where  a   =  density  ratio, 

p  =  standard  sea  level  density 
Th  =  normalized  maximum  thrust, 
HT  =  tropopause  altitude, 

and  h  =  normalized  altitude. 
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Substituting  equations  (A. 57)  thru  (A. 59)  into  equations 
(A. 55)  and  (A. 56),  we  obtain  the  revised  state  equations 


M  =  £—.  Cosa  -  §  sina «§ —  °M2Cn  (A.60) 

a  a  2W  D 

e 

•        gTh   sina        SPoSa     MC      _   gcosy  ,.    6l) 

y  a         M  2W        alR/L  aM  <A.t>i; 


h  =  §^  sin  y    .  (A. 62) 

HL 


The  states  are  Mach  number  M  and  vertical  flight  path 
angle  y.      The  control  is  angle  of  attack  a. 

The  following  inequality  constraints  are  imposed  on 
the  states  and  controls: 


0    <  M  <  MM  (A. 63) 


<*min  1  «  1  «M    •  (A«61<) 


Thrust  (Th)  represents  normalized  maximum  thrust  for 
the  problem  in  Section  VI.   This  is  given  as  a  tabular 
function  of  Mach  number  and  altitude.   The  lift  and  drag 
coefficients  are  given  as  tabular  functions  of  Mach  number 
and  angle  of  attack. 

Empirical  relations  are  used  for  density  ratio  a, 
speed  of  sound  a,  and  maximum  Mach  number  MM  as  functions  of 
altitude.   These  relations  are  given  in  Appendix  D. 
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The  parameters  considered  constant  for  the  problem 
are  the  gravitational  constant  g,  standard  sea  level  density 
p  ,  wing  area  S,  gross  weight  W  ,  and  tropopause  altitude 

V 

JJ .   The  Mathematical  Model  for  the  Missile  Intercept  Problem 

In  this  section  the  mathematical  model  for  the  missile 
intercept  problem  solved  in  Section  V  is  derived.   An  air- 
to-air  missile  launched  from  an  attacking  aircraft  must 
intercept  a  constant-velocity  target.   The  missile  is 
restricted  to  maneuver  in  a  plane.   The  orientation  of  this 
maneuver  plane  in  three  dimensional  space  is  defined  at 
launch  as  the  plane  containing  the  position  of  the  missile 
at  launch,  the  position  of  the  target  at  launch,  and  the 
velocity  vector  of  the  target. 

The  assumptions  applied  to  the  problem  include  those 
presented  in  Section  A.l  plus  the  following: 

a.  the  initial  velocity  vector  of  the  missile  lies  in 
the  maneuver  plane, 

b.  the  attacking  aircraft  is  tracking  the  target  at 
missile  launch  so  that  the  missile's  initial  velocity  points 
at  the  target  at  t  =  0, 

c.  the  target  moves  with  constant  velocity, 

d.  components  of  out  of  plane  forces  perpendicular  to 
the  maneuver  plane  are  ignored. 
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Figure  HO   presents  a  view  of  the  problem  in  the  maneuver 
plane. 


TGT 


Figure  40 
Missile  Coordinate  System 
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Four  axis  systems  are  drawn  in  Figure  kO .      They  are: 

a.  (X',Y')    a  fixed  inertial  axis  system  in  the  maneuver 

plane  with  the  origin  at  the  missile  launch 
point; 

b.  (X  ,Y)     a  Newtonian  reference  system  in  the  maneuver 

plane  with  the  origin  at  the  missile  launch 
point  at  t  =  0;  after  launch  the  origin 
remains  fixed  with  respect  to  the  target 
(it  moves  with  velocity  v„  with  respect  to 
the  X'Y'  system);        l 

c.  (x",y")    a  non-rotating  axis  system  fixed  to  the 

missile  center  of  mass; 

d.  (x',yf)    a  rotating  axis  system  fixed  to  the  missile 

center  of  mass;  the  x'  axis  is  oriented  in 
the  direction  of  the  missile's  velocity 
vector. 

The  systems  X'Y*,  XY,  and  x"y"  are  oriented  in  the  maneuver 
plane  so  that  the  axes  O'X',  OX,  and  ex"  form  the  intersection 
of  the  maneuver  plane  and  a  horizontal  plane.   These  axes  are 
chosen  so  that  the  target's  initial  X,  X',  and  x"  positions 
are  positive.   The  axes  O'Y',  0Y,  and  cy"  are  chosen  so  that 
the  component  of  missile  weight  inthe  maneuver  plane  is  acting 
in  the  negative  Y',  Y,  or  y"  direction.   All  angles  are 
positive  as  they  are  shown  in  Figure  kO   in  the  counter- 
clockwise direction.   The  remaining  notation  is: 

a.  forces; 

N  =  normal  aerodynamic  force  , 
A  =  axial  aerodynamic  force, 
T  =  thrust, 
Th  =  normalized  thrust, 

W  =  component  of  missile  weight  in  the  maneuver 
plane, 

b.  angles; 

a  =  angle  of  attack^ 

6   =  missile  flight  path  angle, 

Y  =  target  flight  path  angle, 
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c.   other  quantities; 

v  =  missile  velocity, 

v™  =  target  velocity, 

M  =  missile  Mach  number, 

MT  =  target  Mach  number, 

W  =  missile  gross  weight, 

C„  =  normal  force  coefficient, 

C.  =  axial  force  coefficient, 

w  =  angular  velocity  of  the  x'y'  system  with 
respect  to  the  x"y"  system, 

m  =  missile  mass 

g  =  gravitational  constant, 

S  =  missile  wing  area 

p  =  air  density, 

n  =  load  factor, 

a  =  speed  of  sound, 

The  equations  of  motion  are  derived  following  the  methods 
used  in  Section  A.l.   Equations  (A.l)  thru  (A. 5)  are 
identical.   The  angular  velocity  of  the  x'y'  system  with 
respect  to  the  non-rotating  frame  x"y"  is  given  by 


w  =  6  e  ,  .  (A. 65) 


The  product 


a)  x   v  =  v  6  e  ,  (A. 66) 


is  formed.   Summing  forces  in  the  x'  and  y'  directions,  we 
obtain  from  equation  (A.l) 
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Tcoscx   -   Nsina   -  Acosa   -  W  sine    =   mv      ,  (A. 67) 

Tsina  +  Ncosa  -  Asina  -  W  cose  =  mv  9   .  (A. 68) 

c 


Axial  and  normal  force  coefficients  are  defined  by  the 
expressions 


A  =  CA  I  pv2S   ,  (A. 69) 


and 


N  =  CN  I  pv2S   .  (A. 70) 


Substituting  equations  (A. 21),  (A. 22),  and  (A. 23)  into 
equations  (A. 69)  and  (A. 70)  and  transforming  the  results 
into  state  variable  format,  the  state  equations  become 

M  =  ^3  Thcosa  -  gp  Mucosa  -  §p  M^sina  -  ^  sine  ,     ( A  .  71 ) 
e  e  e  e 


'  =  ^M  Thsinq  _  ggSa         ggSa       _  ^c  cose.       ,   ?2) 

9   aW     M     2W   il0ASina   2W   U0NCOSa   aW   M  '      ^.^; 
e  e  e  e 


Two  additional  state  equations  are  required  to  impose  end 
conditions  on  the  relative  positions  of  the  missile  and 
target  in  the  optimization  procedure.   They  are 


X  =  aMcos6  -  aMTcosY  (A. 73) 


Y  =   aMsine   -   af^siny      .  (A.71*) 
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The  states  are  missile  Mach  number  M  and  missile  flight 
path  angle  8.  The  control  is  missile  angle  of  attack  a 
Normalized  thrust  is  given  by 


Th  =  1         t  <  tB  (A. 75) 


Th  =  0         t  >  tB 


where  tR  represents  engine  burnout.   The  following  inequality 
constraints  are  imposed: 


"aM  1    a        1  aM  (A'77) 


-nM  1  |  (6M)  <  nM  (A. 78) 


Equation  (A. 78)  represents  a  structural  load  factor  limit. 

The  axial  and  normal  force  coefficients  are  given  as 
tabular  functions  of  Mach  number  and  angle  of  attack. 
Parameters  considered  constant  for  the  problem  are  the 
gravitational  constant  g,  maximum  thrust  ThM,  the  speed  of 
sound  a,  missile  weight  W  ,  air  density  p,  missile  wing  area 
S,  the  Mach  number  of  the  target  M™,  and  the  flight  path 
angle  of  the  target  y. 

In  order  to  properly  define  the  problem,  it  is  necessary 
to  perform  several  manipulations  in  analytic  geometry. 
First,  the  three  dimensional  positions  of  the  missile  and 
target  must  be  specified  at  launch.   Second,  the  velocity 
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vector  of  the  target  and  the  Mach  number  of  the  missile  at 
launch  must  be  specified.   Once  this  is  done  it  is  necessary 
to: 

.a.   identify  the  maneuver  plane, 

b.  identify  the  XY  coordinate  system, 

c.  calculate  the  target  coordinates  in  that  system, 

d.  calculate  the  target  flight  path  angle  y>  the  initial 
missile  flight  path  angle  0(0),  and  the  component  of  the 
missile  weight  acting  in  the  maneuver  plane  W  .   The 
optimization  procedure  can  then  be  commenced. 

To  accomplish  these  calculations  an  initial  coordinate 
system  is  established  in  which  the  problem  can  be  easily 
visualized.   The  origin  is  situated  at  the  missile.   The  OX 
axis  is  positioned  in  the  horizontal  plane.   The  0Y  axis  is 
positioned  in  the  horizontal  plane  such  that  the  target  has 
no  Y' coordinate .   The  0Z_  axis  is  positioned  in  the  vertical 
plane  such  that  a  target  which  has  an  altitude  advantage 
over  the  missile  has  a  positive  Z_  component.   This  coordinate 
system  is  shown  in  figures  4l  and  42.   The  angles  6T  and  &T 
are  defined  as  shown  above.   The  following  relations  may  be 
written: 


a  =  RTi  +  hTk  (A. 79) 

MT  =  mx,i  +  mylj  +  mz,k  (A.80) 


m 


t 


tane^  =  -£'  (A. 81) 

1        mx' 
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Figure  k2 
Initial  Missile  Coordinate  System 
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m  , 
tanB   =  _l_  .  (A. 82) 

m_. 


With  the  problem  defined  in  the  coordinate  system 
shown  in  Figures  4l  and  42  it  is  now  necessary  to  transfer 
the  problem  to  the  coordinate  system  used  in  the  optimization 
procedure.   That  is,  it  is  necessary  to  identify  the  maneuver 
plane  and  the  XY  coordinate  system.   To  this  end,  a  vector 
normal  to  the  maneuver  plane  is 


N  =  a  x  MT  (A. 83) 


=  -hTmy,i  +  (hTmx!-RTmz,  )j  +  Ityryk   .  (A. 84) 

To  establish  the  X  axis  a  vector  is  required  which  is  in 

both  the  maneuver  plane  and  a  horizontal  plane.   Such  a 
vector  is 

X  =  N  x  k  (A. 85) 

=  (hTmxt-RTmz, )i  +  tyn  ,j   .  (A. 86) 

To  establish  the  Y  axis  a  vector  is  required  which  is  in 
the  maneuver  plane  and  perpendicular  to  the  X  axis.   Such 
a  vector  is 
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Y  =  N  x  X 


-RThTmh ,  1  +  RTmy , (hTmx , -RTmz , ) j 

2    2  2 

=[hT  m  ,  +(hTmx,-RTmz,)  ]k   . 


(A. 87) 
(A. 88) 


The  angle  $   between  the  maneuver  plane  and  a  horizontal 
plane  is  required  and  is  given  by 


cos<{> 


N  •  k 


jN|  |k| 


(A. 89) 


RTmy' 


C(hTm   )   +  (hTmx,-RTmz,)"  +  (RTm  ,)  3 


(A. 90) 


The  missile  weight  component  in  the  maneuver  plane  W  may 
be  found  from 


W  =  w  sin<f> 
c    e   y 


(A. 9D 


This   is   shown  graphically   in  Figure    43 


maneuver  plane 


horizontal  plane 

Figure   43 
Missile  Weight   Component 
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The  initial  target  coordinates  are 


xT(o)  = 


PROJxa | 


RT(hTmx,-RTmz, ) 


[(hTmxl-RTmzt)2  +  (hTmy,)2]^ 


(A. 93) 


YT(0) 


±PROJ„a 


(A. 95) 


=  + 


2      2    3    2  2 

-RT  hrpm  ,  -hT  my,  -hT(hTmx,-RTmzt  ) 


The  sign  of  Y„(0)  is  resolved  by 


if  hT  _>  0  ,   YT  is  positive; 


if  h~  <  0  ,   YT  is  negative. 


The  initial  missile  flight  path  angle  9(0)  is 


tane(0)  = 


YT(0) 
J4X0T 


(A. 96) 


by  assumption  b  at  the  beginning  of  this  section.   Before 
proceeding  it  is  necessary  to  insure  that  the  X  and  Y 
vectors  given  in  equations  (A. 86)  and  (A. 88)  have  the 
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correct  sense.   This  may  be  checked  by  observing  the  sign 

of  PROJYa  which  should  be  positive  and  the  sign  of  PROJva 

A—  X  m 

which  should  be  positive  if  hT  >  0  or  negative  if  hT  <  0. 

After  the  senses  of  these  vectors  have  been  checked  and 

altered  as  required,  the  target  flight  path  angle  y   may  be 
calculated  by 


cosy  =     " —  •  (A. 97) 

|M  | |X| 


The  possible  range  of  y   is 


-it  <   y  <   it   .  (A. 98) 


If  cosy  is  positive,  then 


-  \   <  y  <  J   .  (A. 99) 


If  cosy  is  negative,  then 

|  <  Y  <  tt   or   -ir  <_  y  <   -  |  .  (A. 100) 

To  find  which  inequality  applies  in  (A. 100)  the  quantity 

k  =   ~X  "  (A. 101) 

lSj.Mil 

is  formed.   If  k  is  positive,  then 

0  <  y  <  u     •  (A. 102) 
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If  k  is  negative,  then 

-it  <  y  <  0   .  (A. 103) 

This  logic  completes  the  set  up  of  the  problem  in  the 
maneuver  plane. 
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APPENDIX  B 
TABULAR  FUNCTIONS 

In  this  Appendix  the  tabular  functions  used  in  the 
problems  are  presented. 

1.  Three  Dimensional  Plots 

Three  dimensional  plots  of  all  tabular  functions  are 
presented  here.   Figure  44  gives  the  lift  coefficient  Cy  as 
a  function  of  Mach  number  M  and  angle  of  attack  a  for  the 
the  supersonic  fighter  aircraft  used  in  the  aircraft 
problems.   Figure  45  gives  the  drag  coefficient  C~  as  a 
function  of  Mach  number  M  and  angle  of  attack  a  for  the  same 
fighter.   Figure  46  gives  normalized  maximum  thrust  Th  as  a 
function  of  Mach  number  M  and  altitude  H  for  the  supersonic 
aircraft  performing  the  minimum-time  climb  in  the  problem  in 
Section  VI.   Figure  47  gives  the  normal  force  coefficient  CN 
as  a  function  of  Mach  number  M  and  angle  of  attack  a  for  the 
air-to-air  missile  used  in  the  problem  in  Section  V.   Figure 
48  gives  the  axial  force  coefficient  C.  for  this  missile. 

2.  Tables 

Following  each  plot  a  condensed  version  of  the  table  used 
in  the  computation  is  presented. 
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Figure  kH 
CL  =  f(M,a) 
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Figure   45 
CD  =   f(M,o) 
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Figure    46 
Th  =   f (M,H) 


177 


\     \ 


< 

cc 

o 

O 

cc 

O 

HH 

cc 

LO 

«I 

LU 

• 

co 

O 

o 

2TLU 

l-H 

=>Q 

■z. 

ZD 

o 

»- 

</> 

X-< 

o 

cc 

Oh 

o 

LU 

•d-J 

<4- 

o. 

s:< 

• 

z> 

o 

to 

n  ii' 

< 

ora: 

CO 

r£ 

UJLU 

h-i- 

o 

o 

LULU 

o 

LL 

XS 

m 

LU 

<< 

• 

_J 

1— 

O'.CC 

o 

CD 

I/) 

<< 

< 

3 

O-Q, 

K 

CC 

X 

_J 

h- 

< 

o 

l-_J 

a 

3T 

2T<f 

(M 

3 

Oo 

• 

2u 

rvjt-i 

o 

H-« 

»-<l- 

X 

c£i£ 

< 

OLU 

s: 

x> 

o 
o 

0<M>OO^r»-ir>or->0'-fc\JNj-  oo>«i-'-<or<"i<\i>i->i-oOf>'-r\j 
O  cor^-J-0N<M<MyN'^fM'-HO^-ir\'^(NicorOrorOsO<7^— <mCT*m 

O         cr rr> o --* -4" r^oiAo^>i-o>o ro ^cj^r^  .o>n  -t-coiMoj.- ioo 

00  cOM^^^^-tmMfviN^^^OOOOOOOOOOO 

o        ooooooooouooooooooooooooo 


O  >-«  rsi  o  O  co  *-*  -o  *-*  •$■  <?  vO  r»-  -O  oo  oo  cm  co  o>  -J-  in  co  ~0  in  — < 

o        f\j  J- ^oor- cor*- ^-j- -^vj-o --<  co  o^tn '-!•-•■-•  loco— <^o  m 

O         r-r~cootNJincr -j- corn  ^mroo  <x>  r-»voir> -4- mrgog^Ho  o 

r-        cor-.ooin<}-i'ococMoj^-<-Hi-i-Hooooooooooo 

•    •  •#§••••••••••••••••••••• 

o   oooooooooooooocoooooooooo 


o^r^inrnco^ocorooococoino^oor--ino>inoNr~cMCMOv 

o        ^-iin^-r*- vj-r*-r-ror»-cocMncOso^^'c^ooNj-r>-o  j-o^^- 

o        minso:0'-<^"0ocoi^i>J"om.NjO"3^~u"m"i^-rt.(\j%i--'00 

vO        cor-.ommvt.firoiNJcMr-t^r-i^ooooooo  ^ooo 

•    #•••••••••••••••«•••••••• 

o   ooooooooooooooooooooooooo 


Ov0ooroco^c7^cr,'-,,^inr0v0f~-'^ro^r<roooc\jro0f>oco 
cm  co  r\j  -^  cm  t^  cc  ir  \  ^  co  in  cm  in  ,-o  n  cm  co  ■>  o  co  r~  O  '•o  o^  -J" 
co  com  i*»  o  «<ir—  cm  r-  cm  co  m  cm  o  cor-  in  -j-  »?■  m  -m  ai— <  o  o 

(UN>0,^in<TI,lffl'MN^HrHHOOOOOOOOOOO 

ooooooooooooooooooooooooo 


o  r- m -0  lo  ^-i  i  o -^  r~- o  cm  O  "i  (\i  <f  rW~- rn  o^  ^ -J3  O  rn  cri  nJ- 
cm  cm  ^t- vO  (^  m  r- rg -o  cm  co  in  cm  o  *  p>-  a  v  j- m  co  cm  cm  ^  o  o 
a>  r-  >o  m  xj-  <r  co  ro  <\i  "M  h  •-<  ■-«  — <  o  o  o  o  o  o  o  o  o  o  o 

OOOOOOOOOOOOQOOOOOOOO'-'OOO 


O  o  v0  cm  in  cm  O  «■>  o*  co  vO  cm  co  o>  4-  vj-  ct*  --i  o  CMn  r-  in  co  r— 
oomo^vOcors-cor>-ococMOroo-ODcofMvOC?,<,"ico^" 

r-l— ( CO  in  CO  CM  -O  — tO^(CC-tCMOro>-ir,  -^-  T\  rO  CM .- <  — <  O  O 

CONvOinJ-'trtl^M.M^-li-iHOOOoOOOOOOO 

ooooooooooooooooooooooooo 


oco.^cocoinN.in^<>.^oco^r^co<f  r-nj-oco  ^•-J-^-r-■- 
cocMO^t•lnoJln^J••--'lnv>r»'-■^ocM:^■or--cocM>oo'Mnco>i■ 

O  ^-<  CO  in  oo  CM -0 -^ -O  ^H  r— -4"  CM  o  CO  o  m  ^J- CO  .-o  CM  •-<>-<  3  o 

(rNOn-t-l-MtfiNM^-ti-i-ioOOOOOOOOOO 

ooooooooooooooooooooooooo 


o  oo  — « co  co  Ln  r»  m  -o  o>  <-\  o  co  —•  r«-  co  ■$■  r»  o  o  co  o  -t  r»  r- 
cO'MOvj-Lnrsimvj-— <inJsr---<ocMj>vOr---~3cM^aMOco^ 
O^co^cocmO'HO— tr-sj-CMOco-onvTcoco^J— moo 

MNO'n-J-vfCifflMnjHHHHOOOOOOOOOOO 
OOOOOOOOOOOOOOOOOOOOOOOOO 


o  o  -o  rg  m  cm  o^  co  o  co  o  cm  n  cr>  -j-  <f  o  ^  o  o  m  r-  m  co  r~ 

ooMno>cr>-oor-cor-ococMO>-ooocococMOO>((ino^- 

r-tHfvricoMO'-'~o-<w-}'iMO'or»  n  >r  ■"'I  co  cm  r-i  r-i  o  o 

O  OOr-OLO-J-^rrOCO'MCN.-l— <  -«— 100000000000 

•  •##•••#•••••#•••••••••••• 

o   ooooooooooooooooooooooooo 


ooooooooooooooooooooooooo 

ooooooooooooooooooooooooo 
oooooooooooooooooooooooo 

OOOOOOOOOOOOOOOOO  3000000 

inoinowomoinomoinoioouVT^ouiouio 
.-i  •"-!  cm  cm  ro  co  <f-  .4-  m  in  -0  -o  r-  Is-  oo  co  0^  cr>  o  O  >-t  ^  im 


178 


\ 


oir\vca>^"Cr*<x>ir\Lnr-o  -orno^-NO  -4-c0v0r-irMr~of-m 

O         Ococoo<iHM^rooa)^ooK01oco-},0'HHc^r<ir>- 

O         <Or-tooO  o  JD^-CT>'-HrncorncMri<MOno^--ol-n^"rrirM'-<0 

h-        rjrtOiooNJDm^^m^MHH-HHOoooooooo 

•    ■•••••••••■••••••■••••••• 

rW  rHr-toOOOOOOOOOOOOoOOOOOOOOO 


o  o^om  — 4~o^ro-j-Lna^r^(\io^'Ln  vj-r-icocjvooomr*- 
O  »h  r^  vr  en  ro  co  in  f^  ^  cm  no  cm  coin  cm  o  cor- in  ^rco  cm  cm  r^o 
■O        (^o^wsom^rfimrgiMHHMHOoooooooo 


• 


< 

a: 

o 

0 

a 

0 

•—1 

Ql 

^ 

< 

UJ 

• 

CD 

r-l 

o 

SUJ 

•a 

0 

•H 
■P 

c 
o 
o 

►— 1 

DQ 

Z 

ZD 

a 

h- 

CO 

I- 

O 

cC 

Ot- 

O 

UJ 

<t_l 

f0 

Q. 

^.< 

• 

r> 

r-l 

00 

H   II 

1 — 

< 

ccai 

m 

ujuj 

iH 

a£ 

\~\- 

O 

a 

UllU 

O 

u. 

2:5: 

<M 

UJ 

<< 

i 

_J 

H- 

arnc 

r-l 

03 

co 

<< 

< 

Q.Q. 

X 

_J 

h- 

< 

O 

h-_l 

O 

s: 

Z< 

.-I 

=) 

OO 

t 

s: 

Mr-. 

r-l 

»— « 

->h- 

X 

a:cy. 

< 

OLU 

z. 

x> 

O 
O 
O 

• 

1-1  r-H  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 


o  oooOr-ion)~t'X)-HffKMin-to>OH'Dmr»(Da)>(\jo 

o  -ocNjr-too— i'Mmr-r-nnr-H|— Nj-r-toco^in-JTcvM-Hr-io 

in  .-ioo'X>r-~H3mNt<<v,"><MCM.-ir-ir-irwooooooooo 

•  •■••••••••••••••••■•••••• 

—*  rHr-iOOoOOOOOOOOOOOOOOOOOOOO 


o  <^  rH  o  rg  r-^  ro  ro  o  in  >j- cr  c\i  00  r~  cr  (\i  r-t  tn  <\i  1a  eg  «o  ^  ir> 
•4"oovj-Ot^i-or^->ofriox)"fi'X)3o  -j-  oco-oroiTi  vOi~-coc\i  o 
r-icor— ^or^a>oc<">vou>>ro  vjfir-iuvr--  .om-r  mcM-H^o 
r-4^aor-  oinif,sj-m,MiMi\jHHHOooooooooo 

r-l  -H00000O0O0OOO0OOOO00OOO03 


O  ( T>  sO  O  vO  O  rO  —1  nC  r-  •£>  no  v0  r<^  —I  O  i-i  in  -J"  4"  O  —*  CT>  -O  CM 

00-  'J^to  >oc^r^c>or^coinr-iroo''<vnl-o— irotn>or-r-is0 
r^^-rn.~o  ^<3corH«^a3coo^vOc<"i--ic/»r--soin>}-rn<MrHrHO 

ooxaNjjinvj-^ffiiNj'MHHHHaooooooooo 

r-<        r-noooooooooooooooooooooooo 


o  m  cm  o  .4-  a*  r-i  -J-  o  o  Is-  m  in  m  o  -o  rn  cm  in  00  r~  — t  cm  <m  o 
o>l-fO'--<cM'--<cjM<*i'-n-o(^a5m<.,Os0^f^<-<o^r-<ro-nr»'-io 

Mr-(OrH(M>J-J)Oii1|^fMC0inMO3)N  O-t-tflMMMO 
O(M0r^^>u"\<f  ^iWixj^hhhOoOOQOOOOO 

—1000000000000000000000000 


Oror^r-oocOr^cMn— i<MfM>^inr-CMLn-l-orooinNN 
~omrooincoo>NO'-ima>0'OrM— iino>'W~ocM<)-^o.n 
corx-r--roOr-i<t-co(MvO'-<cx)<fcMOCco!n^-^-noj— mo 
cjmo>"-  oinin-j-romcMfM— <— t^Hr-io  0000  00000 

r-l     OOOOOOOOOOOOOqOOOOOOOOOOO 


o  in  o^  o  co  in  co  00  00  c  •")  4- cm  Ln  r\i  t>  cm  cm  co  o  co  rr»  rn  o  >4- in 

r^ooocornu^-j-r-ir^rnin^  jx;oror^vO  oco-*m-i)Oi^ 
in-tin-or~oi'^r~r-iu"\'-<r--'^--io^ooNon-4-rorocMr-i--4o 

OiCOS<Jif>invf^rflNiNr-tHHOOOOOCOOOOO 
rH  OOOOOOOOOOOOOQOOOOOOOOOOO 


o  cm  o  cj>  -4-  cm  co  (/>  r-<  >j-  0  —1  in  r-<  —( {f\  o\  cr>  -j-  m  cm  in  -j-  o  -j- 

o        .^-O'.n^-'-o  oOrH^^j^^in^^-j-i-siT^ON-nOb-* 
O         cm  cm  cm  r<">  m  co  cm  ^J  0^0  Oct  r-l  cr>co  -om<r  rnr<ii\MrHo 

CT<  0>COr~-vOm  sfvJ-cri-T)(MCMr-i— l-HOOOOOOOOOOO 

O  OOOOOOOOOOOOOOOOOOOOOOOOO 


OOOOOOOOOOOOOOOOOOOOOOOOO 

OOOOOOOOOOOOOOOOOOOOOOOOO 
OOOOOOOOOOOOOO  OOOOOOO  CO  00 

OOOOOOOOOOOOOOOOOOOOOOOO 

in  o  m  o  1  n  o  m  o  in  o  in  0  m  o  u">  o  u  \  0  m  o  m  o  m  o 

t-i  rH  cm  im  o  *  ro  -4-  >j-  m  in  o  -o  i*-  r*-  co  co  cr>  <j>  o  o  .-h  —i  <m 


179 


\ 


o 

o 

•4- 

• 

rg 

V- 

ll 

< 

C£ 

o 

o 

ex 

o 

»— 4 

q: 

ro 

< 

LU 

• 

CD 

cm 

o 

SZJJ 

■a 

0) 

c 

•H 

r— t 

oo 

z 

zz> 

o 

t- 

to 

IM 

o 

Oi 

<_>!- 

o 

1X1 

<t-l 

r\J 

c 
o 
o 

Q. 

i;<; 

• 

3 

C\J 

00 

ii  ii 

* — 

< 

oca: 

CO 

rH 

1* 

LULU 
HI- 

o 

o 

LULU 

o 

u. 

sis: 

-H 

LU 

<< 

• 

— 1 

i- 

or  a: 

rj 

co 

OO 

<< 

< 

3 

Q.Q. 

\- 

oc 

X 

_l 

h- 

< 

o 

H-* 

o 

X 

Z< 

o 

Z> 

ou 

• 

S 

rsi.-H 

rg 

l-H 

•-•1— 

X 

a-vC 

< 

OLU 

s 

I> 

o 
o 

o  o  r~  vj- -t  h- oo  oo  r\j  r^- rn  o  cm  ^-i  o  r- o  ia  (\j  >j- o  r~  >n  ia  ct> 

Oc&--tOi0t7>C0'^r--C>r-<rri00r0v£i00O— «~slO^  -O»-i"OC0.> 

—i  t-i  -j-  oo  c<">  .?>  >•  vu  m  s\  co  ■— •  m  i-i  r~  -T  r\i  o  oo  -o  ^"»  -r  <ni  -h  o 


•OOOOOCOOOOOOOOOOOOOO 


O'-KMf-'TNaoo  <5'v-<"<^aD  t  — <v0>ro  ror-noo^rocom 
comcarnco  o^T-rrn-r -o^^Ovj-Tr-<i/'r^Nr>inrn<N'^o 

^H^r-i^OOOOOOOOOOOOOOOOOOOOO 


o  r-  cr*  co  co  —<  o  o  n  —<  — i  f*>  •-*  ^o  -o  m  <\i  r g  ^-  r-  »j-  ro  cm  o  —« 
o  fiH.<ia>oinir,MNO^-NLiH  oOt'OmrOi-icoor*-^ 
h-  lt>  r<v -o  -J-  f\i  --*  *-i  —(  cm  m  oo  i  n  o*  -o  ro  —» cr*  r-  -o  m  rocvj  *-*  o 
inmrvJU^^I^vOin-TrOMM-i-ii-irHOOOOOOOO 

r-tr-ir-4— iooooooooooooooooooooo 


o  v0  o  in  in  <m  i-i  m  Is-  ro  in  n  in  <T>  -J-  -4-  m  r-i  r-i  o>  h  -o  o  ro  p~ 
Or-if^mvoo^'or^-coro  4*^*>o  .o«i-oino<\jOL>vom>oco 
Or^NflO  »C0f0c0Oc0'<-<M00mr<,iOf3r-«0'4-'*>'VJ— <o 

'H'H-ir-ioOOOOOOOOOOOOOOOOOOOO 


OC>0'0(NjOir\HvOrOOOMJ3ScoM\IO(NOOOOm 

in  rn  v0  o  r- in  r<"K\j  r*- in  o 'M  o  cd  r~  4- o  x\  qn  33  f- m  •  j- in  co 
rni*-<*\jO  0ir\LOv0vOixi(\i  o<-^r-  4-^JOco-om*tm<\j.-<o 
<t(\ir-icr>  BNvOin^ntn^iM'-f^HHOooooooo 

r-H^lr^OoOOOOOOOoOOo^OOOoOOOO 


ooo^'Momovi oOhvON  or^inonor-tcoomoo^o 

mOX>r*O<'rlv0X>c0-J>s0^f^'-<,-,J'>v0  r-UDinu^rOrOvJ-00 

r^-cM^--^- r>c^jcNjro4-vOO'n'jrs-<r'-^ J^  •.oom^rrKvj— <o 

(•OCMOCJncOI^- -Ou^J-i^i-oN.MrHi-i-cDOO  DOOOOQ 
•-H  r-<-Hn-lOOOOOOOOOoOOOOoOOOOOOO 


o»-i^— 'Ocr>c0'JN'-<'-<'^O'3v'— '-0  in  <j*^j"0  >!■—<— 'Oi"no 
o         •-i>oroocru>o^<Hc\jiniJ>-J,o^^Jf<'i'^u<r,-om<rr>jtMi-<u 

00  WHOl>N'Uin'n'T(flM!\HHr<HOOOOOOOOO 

,-i  HpIHOOOOOOOOOOOOOOOOOOOOOO 


ooooooooooooooooooooooooo 

OOOOOOOOOOOOOOOOOOOOOOOOO 

oooooooooooooooooooooooo 

OOOOOOOOOOO  JOO'3  0  000  000  00 

inomoino^ouio^oinoiAOiriOinouyo^O 
»-(■-( i%J(Mr0r>ivi-  ^■>Ain'0>Ol^r»vOa3^Lr'00-ir-i(M 


180 


Figure    A  7 
CN  =   f(M,a) 
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Figure    48 
CA  =   f(M,a) 
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APPENDIX  C 
INTERPOLATION 

In  the  optimization  problems  solved  herein  the  aero- 
dynamic data  is  given  in  tabular  form.   The  dependent 
variable  D  is  given  as  a  function  of  two  independent 
variables  M  and  a  in  all  cases.   Excerpts  from  the  tables 

used  in  computation  are  presented  in  Appendix  B. 

2      2 

3D    3D    3D    3D 
For  a  given  M  and  a  quantities  D,  -ttt  ,  ~—  ,  — »  ,  — x   , 

3M    3a 

2 
3  D 

and  .„-    are  required  by  the  optimization  algorithm. 

Parabolic  interpolation  is  used  to  obtain  these  quantities. 
In  this  Appendix  parabolic  interpolation  for  two  independent 
variables  is  derived. 

1.   Parabolic  Interpolation  in  the  Plane 

To  apply  parabolic  interpolation  to  a  tabular  function 

of  two  independent  variables  the  nearest  point  in  the 

tables  to  the  given  point  (M,a)  must  first  be  found.   It  is 

assumed  that  the  tabular  data  is  given  at  constant  intervals 

AM  and  Aa  in  the  independent  variables.   The  nearest  point 

given  in  the  tables  and  the  surrounding  eight  points  are 

required  in  the  interpolation  and  are  shown  in  Figure  49. 

The  parameters  6  and  <j>  locate  the  point  (M,cx)  from  the 

nearest  tabular  point  (M  ,a  ).   If  (M  ,a  )  is  the  nearest 

*  s*  s         s'  s 


point  then 


|  <  4>  <  |  (c.i) 
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(Ms-rV 


(Ms.ra 


) 


<Ms'as+ 


Aa 


JL 


)         (Ms,os) 


<MS-Wl>      <Ms'V! 


<{)Aa 
t  (M,a) 

6  AM 


AM" 


) 


(Ms+rVi} 


<Ms+ras) 


<\n-Vi> 


M 


Figure    4 9 


and 


*!•!* 


(C.2) 


The   inequalities    (C.l)    and    (C.2)    hold  unless   the   nearest 
point    (M    ,a    )    is   on  a  border   of  the   table.      In  this   case 

5         S 

(M  , a  )  is  chosen  one  point  in  from  the  border.   The 
s   s 

parameters  9  and  0  then  satisfy 


-  1  <  <J)  <  1 


(C.3) 


and 


-  1  <  6  <  1 


(C.I) 
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Writing  Taylor  series  expansions  including  up  to  second 

order  terms  for  the  eight  points  surrounding  (M  ,a  ) ,  we 

s   s 

obtain 


D(Ms+AM,as+Aa)  =  D(M8+1,as+1)  ■  D(Ms,as) 


(C.5) 


+  ^W 


+  Aa  ^~ 
9a 


+  (AM)2  32D    +  (Aa)2  32p 


2    3M2 


2    3a2 


+  AMAa 


32D 


D(M  -AM,a  -Aa)  =  D(M   , ,a   _)  *  D(M  ,a  ) 
s    '  s  s-1'  s-1       s'  s 


(C.6) 


-  AM 


3D 

3M 


-  Aa  ;r— 

3a 


+  (AM)2  32D 
s      2    3M2 


(Aa)2  32D 


3a' 


+  AMAa 


32D 

3M3a 


D(M  -AM, a  +Aa)  =  D(M   . ,a  .. )  =  D(M  ,a  ) 
s    *  s  s-15  s+1       s*  s 


(C.7) 


-  AM 


3D 

3M 


+  Aa|° 
3a 


+  (AM)2  32D    +  (Aa)2  32D 


2    3M2 


3a2 


-  AMAa 


32D 


D(Ms+AM,as-Aa)  =  D(Mg+1,as_1)  *   D(Ms,ag) 


(C.8) 


+  AM 


3D 

"3~M 


-  Aa 


3D 
3a 


+  (AM)2  32D    +  (Aa)2  32p 


2    3M2 


3a2 


-  AMAa 


32D 
3M3a 


3D 


D(M  +AM,a  )  =  D(M  ,.  , a  )  =  D(M  ,a  )  +  AM  ^ 
s    '  s       s+1'  s       s*  s       3M 


+  (AM)2  32D 
s      2    3M2 


(C.9) 
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D(M  -AM,a  ) 
s     s 


D(Ms_1>a5) 


D(Ms,as) 


AM  & 
n   3M 


.  (AM)2  32D 

2      2 
s      d        3M 


(CIO) 


D(M  .a  +Aa) 

S    S 


=  D(Ms,as+1)  " 


^V0^  +  Aa  f£ 


(Aa)2  32D 

2    a  2 
s         3a 


(C.ll) 


D(Ms,as-Aa) 


D(Ms,as_1) 


D(Ms>as) 


a   3D 

Aa  s — 

3a 


+  (A")2  32D 

2    »  2 
s         3a 


(C.12) 


Subtracting  equation  (CIO)  from  equation  (C.9),  we  obtain 


D(Ms+l*as)  "  D(Ms-l>as)  '  2AM  I 


3D 

M 


or 


3D 

3M 


D(Ms+l'as)  -  D(Ms-l'as) 


2AM 


(C.13) 


Subtracting  equation  (C12)  from  equation  (C.ll),  we  obtain 


3D 


D<*B«W  "  D(MS'°B-1)  =  2Att  IS 


or 


3D 
3a 


D<Ms'°s+l)  "  D(Ms'as-l) 


2Aa 


(C.1H) 
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Adding  equation  (CIO)  to  equation  (C.9),  we  obtain 


D(Ms+l>as)  +  D<Vl-«.>  =  2D(Ms.«s)  +  (AM)2  "~3 


3M 


or 


92D   «   D(Ms  +  l'«s)  •  2D<Ms'as>  +  ^"s-l^s5 


8M 


(AM)' 


(C15) 


Adding  equation  (C.12)  to  equation  (C.ll),  we  obtain 


D(MB'a8  +  l)  +  ^S^S-l5  "  2D<MS'°S>  +  (Aa)2  H 

3a 


or 


92D 


8a 


D(Msiqs+1)  -  2D(Ms>as)  +  D^a^) 
(Aa)2 


(C.16) 


Subtracting  equation  (C.7)  plus  equation  (C.8)  from  equation 
(C.5)  plus  equation  (C.6),  we  obtain 

D<*W<W  +  D«Vl«°8-l)  -  D(MS-1'°S+1)  -  D(MS+1'°S-1)  Z 


IJAMAa 


92d 


9M9a 


(C.17) 


or 


3M8a 


JTAMAa 
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A  Taylor  series  is  written  for  the  point  (M,a).   Using 
equations  (C.13),  C.ll),  (C.15),  (C.16),  and  (C.17),  we  may 
reduce  this  series  to 


D(Ms+eAM,as+<f>Aa)  =  D(M,a)  =  D(Ms>as)  +  |  [D(Ms+1,as)  -  DO^,^)] 

+  I  ^^s'W-^s'Vl"  +  V  ^(MS+l»0s)-5DCMS'OS)4D(M8-rO8):i 

+  V  tD(Ms,as+1)+D(Ms)as_1)-2D(Ms,as)] 


'+$  ^s+i'Vi^Vi'Vi^s-i'Vi^s+i'Vi"    ■ 


Rearranging,    we  have 


D(M,cO   =  *$■  D(M8_1,oB_1)   +  iSA^n.  D(Ms,og_1) 
"  ¥  D<MS+l'«s-l'    +  ^T^  D»W».> 


+    (l-e2-*2)    D(M.,o.) 
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Equation  (C.19)  is  the  expression  used  to  interpolate  for 
the  value  of  D  in  terms  of  the  surrounding  none  tabular 
points.   To  obtain  expressions  for  the  required  partial 
derivatives,  observe  that 


M  =  M   +  6AM  (C.20) 

s 


and 


a  =  a   +  <J>Aa   .  (C.21) 

s 


Therefore, 


d9   _1_ 

dM  "  AM 


and 


(C.22) 


(C.23) 


d£  _  _1_ 
da   Aa 


The  chain  rule  for  partial  derivatives  yields 


|S  «.♦•«.•»+♦*«>  -  |§  (M..)  -  Sfig^Sl  ||  .        (C.2M 


Taking  the  partial  derivative  of  equation  (C.19)  with 
respect  to  6,  we  obtain 
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§  CM,a)  -  ^  [  f  DCM^.a^)  -  f  D(M8+1,  a^)  ♦  ^  D(Ms+1,as) 


+  I  D<W<W  "  I  D(Ms-l'<W  +  HT  D(rVl»aS)        (C'25) 


-  29  D(Mg,as)]       . 


Using   similar  procedures,    we  may  derive   the   remaining 
expressions, 

I  (M,a)  -  i  C  f  DCMj ^  ♦  ^  DCM^a^)  -  £  DO^.a^) 

+  !  ^s+l'^W  +  2JT-  ^s'Vl5  "  f  D(Ms-l'<W        (C  ■ 26) 


-  2<|>  D(Ms,as)] 


32D 


(M,a)  =  — ^  [D(Ms+1,as)  -  2D(Ms,as)  +  D^^,^)]  (C.27) 


MT  (AM) 


3: 


2 

£§  (M,a)  =  -^-2  CDCMg.a^)  -  2D(Ms,«s)  +  D(Ms,as+1)]  (C.28) 

3a  (Aa) 


(c.29) 

-  D(Ms-rVi)] 
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APPENDIX  D 
EMPIRICAL  RELATIONS 

In  the  problem  treated  in  Section  VI  empirical  relations 
are  used  for 

a  =  f(H),  (D.l) 

a  =  f(H),  (D.2) 


and 


MM  =  f(H)  (D.3) 


where  the  parameters  are  air  density  ratio  (a),  speed  of 
sound  (a),  maximum  Mach  number  (MM)>  and  altitude  (H) .   The 
air  density  ratio  is  defined  as 

a  -  -fi-  (D.H) 

po 

where  a  is  sea' level  standard  day  density.   This  Appendix 
presents  these  empirical  relations  and  compares  the  values 
obtained  from  these  relations  with  standard  atmospheric 
conditions . 

1.   Air  Density  Ratio 

The  empirical  relation  used  for  air  density  ratio  is 


-c,h      -c0h 
a  =  e  x     +  c  h  *  (D.5) 
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where 

h  =  jf-  (D.6) 

HL 

and 

HL  =  36,089  ft.  (D.7) 

c1  =  1.5^100  (D.8) 

c2  =  1.80445  (D.9) 

c      =  0.4130  (D.10) 


Figure  50  is  a  plot  of  the  values  of  a  obtained  from 
equation  (D.5)  compared  to  those  obtained  from  standard 
atmospheric  tables. 

2.   Speed  of  Sound 

The  empirical  relation  used  for  the  speed  of  sound  is 


a  =  ao(l-c?h)     ,   h  <  1  (D.ll) 

=  971  ft. /sec.  ,   h  >  1  (D.12) 


where  the  parameter  a   is  the  speed  of  sound  at  sea  level 

o 

on  a  standard  day;  that  is 


and 


a  =  1116.89  ft. /sec.  (D.13) 

o 


c?  =  0.1331  (D.14) 


Figure  51  is  a  plot  of  a  vs.  H.   The  expressions  (D.ll)  and 
(D.12)  duplicate  exactly  the  values  obtained  from  standard 
atmospheric  tables. 
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Figure  50 
Air  density  vs.  altitude 
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Figure    51 
Speed  of   sound  vs.    altitude 
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3.   Maximum  Mach  Number 

The  empirical  relation  used  for  the  maximum  Mach 
number  of  the  aircraft  for  the  problem  solved  in  Section  VI 
is 


MM  =  2.1  -  1.1 


-2.l»h 


(D.15) 


Figure  52  is  a  plot  of  equation  (E.9)  along  with  the  actual 
restrictions  of  the  aircraft  under  consideration. 
2.2 


100 


altitude  (ft.  x  NT) 
Figure    52 
Placard  Mach  number  vs.    altitude 
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APPENDIX  E 
A  CONVEXITY  THEOREM 

In  this  Appendix  the  following  theorem  on  convexity  is 
proved. 

Theorem  1;   If  f(x)  is  convex  on  R  where  x  e Rn,  and  f  >  0, 

K  n 

then  f  (x)  is  convex  on  R  where  K  is  any  positive  integer. 

This  theorem  Is  proved  by  mathematical  induction.   First, 

the  following  theorem  is  proved. 

Theorem  2;   If  f(x)  is  convex  on  R  where  x  eR  ,  and  f  >_  0, 

then  f  (x)  is  convex  on  R  . 

Proof:   The  function  f(x)  is  convex  if 


f[Xx2  +  (l-A)^]  <  Xf(x2)  +  (1-X)f(x1)        (E.l) 


for  all  x,,  x~  and  X  e  [0,1],   Squaring  both  sides  of 
Inequality  (E.l)  we  have 


f2[Xx2  +  (1-X)x1]  <  [Xf(x2)+  (1-X)f(x1)]2  .    (E.2) 


The  sense  of  the  inequality  (E.2)  is  retained  as  f  _>  0.   To 

p 
prove  that  f  (x)  is  convex,  it  must  be  shown  that 


f2[Xx2  +  (1-X)^]  <  Xf2(x2)  +  (1-X)f2(x1)  .    (E.3) 


Observing  inequality  (E.2),  (E.3)  is  seen  to  be  a  true 
inequality  if  it  can  be  shown  that 
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[Xf(x2)  +  (1-X)f(x1)]2  <  Xf2(x2)  +  (1-X)f2(x1)  .   (E.H) 


To  show  that  inequality  (E.4)  is  true,  we  proceed  as  follows 
Since  X  e  [0,1],  it  is  true  that 


X[f(x2)  -  f(xx)]2  <  Ef(x2)-  f(x1)]2   .  (E.5) 


Expanding  inequality  (E.5),  we  have 


X[f(x2  -  f(x1)]2  <  f2(x2)  -  2f(x2)f(x1)  +  f2(xx)  .(E.6) 


Rearranging  inequality  (E.6),  we  have 


X[f(x2)  -  f(x1)]2  +  2f(x2)f(x1)  -  f2(xx)  <  f2(x2).(E.7) 


p 

Subtracting  f  (x, )  from  both  sides,  we  obtain 


X[f(x2)  -  f(x1)]2  +  2f(x1)[f(x2)  -  f(x1)]  <  f2(x2)  -  f2(xx)  . 

(E.8) 


Multiplying  inequality  (E.8)  by  X  ,  we  have 


X2f2(x2)  -  2X2f(x2)f(x1)  +  X2f2(xx)  +  2Xf(xx)f(x2)       (E.9) 
-  2Xf2(x, )  <  Xf2(x0)  -  Xf2(x, )   . 
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2 
Adding  f  (x,)  to  both  sides  of  inequality  (E.9),  we  obtain 


X2f2(x2)  -  2X2f(x2)f(x1)  +  2Xf(x1)f(x2)  +  r2(x±) 


(E.10) 


-2Xf2(x1)  +  X2f2(x1)  <  Xf2(x2)  -  Xf2(x1)  +  f2(x1)  . 


Simplifying  Inequality  (E.10),  we  obtain 


[Xf(x2)  +  (1-X)f(x1)]2  <  Xf2(x2)+  (1-X)f2(x1)  .    (E.ll) 


This  is  the  inequality  we  set  out  to  show.   The  theorem  is 
proved. 

Second,  the  following  theorem  is  proved. 

Theorem  3:   If  f  (x)  is  convex  on  Rn  where  xeRn,  and  f  >_  0, 

then  f    (x)  is  convex  on  R  . 

Proof:   It  has  already  been  shown  that  if  f(x)  is  convex  and 

2 

f  >  0,  then  f  (x)  is  convex.   It  is  now  assumed  that 


fK[Xx2  +  (1-X)x1l  <  XfK(x2)  +  (l-X)fK(x1)   .       (E.12) 


Multiplying  both  sides  of  inequality  (E.12)  by  f[Xx2  +  (1-X)x1] 
a  positive  quantity,  we  obtain 


fK+1[Xx2  +  (1-X)^]  <  [XfK(x2)+(l-X)fK(x1)]{f[Xx2+(l-X)x1]} 

(E.13) 
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Substituting  the  expression  (E.l)  into  inequality  (E.13), 
we  have 


fK+1[Xx2+(l-X)x1]  <  [XfK(x2)+(l-X)fK(x1)][Xf(x2)+(l-X)f(x1)] 


(E.lH) 


To  prove  that  f    (x)  is  convex  it  must  be  shown  that 


fK+1[X$2+(1"X)?l]  -   XfK+1(^2)  +  d-^)fK+1(x1)  •    (E.15) 


Observing  inequality  (E.14),  (E.15)  is  seen  to  be  a  true 
inequality  if  it  can  be  shown  that 

[XfK(x2)+(l-X)fK(x1)][Xf(x2)+(l-X)f(x1)]  <  Xf**1^)  +  (l-X)^1^)  . 

(E.16) 

To  show  that  inequality  (E.16)  is  true  we  proceed  as  follows 
The  expression 


[fK(x2)  -  fK(x1)][f(x2)  -  f(xx)]  (E.17) 


is  always  a  positive  number  because  the  signs  of  the  expres- 
sions in  parentheses  in  expression  (E.17)  must  be  the  same. 
The  following  inequality  holds 


X[fK(x2)-fK(x1)][f(x2)-f(x1)]  <  [fK(x2)-fK(x1)][f(x2)-f(x1)] 

(E.18) 
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Expanding  and  multiplying  by  X,  we  obtain 


X2[fK+1(x2)-fK(x1)f(x2)-f(x1)fK(x2)+fK+1(x1)]  < 


X[fK+1(x2)-fK(x1)f(x2)-f(x1)fK(x2)+fK+1(x1)] 


(E.19) 


K+l 
Rearranging  inequality  (E.19)  and  subtracting  2Xf    (x,) 

from  both  sides,  we  obtain 


X2[fK+1(x2)-fK(x1)f(x2)-f(x1)fK(x2)+fK+1(x1)] 

(E.20) 
+  X[fK(x1)f(x2)+f(x1)fK(x2)-2fK+1(x1)]  <  X[fK+1(x2)-fK+1(X;L)] 


Further  expansion  and  rearranging  yields 


X2fK+1(x2)+XfK(x1)f(x2)-X2fK(x1)f(x2)+Xf(x1)fK(x2)-X2f(x1)fK(x2) 


-  2XfK+1(x1)+X2fK+1(x1)  <  XfK+1(x2)-XfK+1(Xl) 


(E.21) 


K+l 
Adding  f    (x,)  to  both  sides  and  rearranging  further,  we 

have 


X2fK+1(x2)+X(l-X)fK(x1)f(x2)+X(l-X)f(x1)fK(x2)+(l-X)2fK+1(x1) 

<  XfK+1(x2)+(l-X)fK+1(Xl)  (E'22) 

or 

[XfK(x2)+(l-X)fK(x1)][Xf(x2)+(l-X)f(x1)]  <  XfK+1(x2)+(l-X)fK+1(x1)   . 

(E.23) 
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This  is  the  inequality  we  set  out  to  show  in  (E.16).   The 
theorem  is  proved. 

It  has  been  shown  that  if  f (x)  is  convex  and  f  <_  0  then 

2  3 

f  (x)  is  convex.   By  Theorem  3,  f  (x)  is  convex.   By 

I] 

Theorem  3  again,  f  (x)  is  convex.   This  reasoning  can  be 

followed  for  all  powers  K  where  K  is  a  positive  integer. 
The  basic  theorem  is  established. 
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Raphson  method  to  minimize  the  epsilon  functional.   The  full 
Newton-Raphson  method  exhibits  terminal  convergence  characteristics 
superior  to  the  "modified"  method,  whereas  the  "modified"  method 
is  generally  superior  in  the  initial  stages  of  a  problem.   An 
algorithm  is  developed  which  uses  both  techniques  in  a 
complementary  way. 

A  new  penalty  functional  which  has  desirable  theoretical 
Droperties  and  exhibits  excellent  computational  behavior  is 
introduced  to  treat  state  and  control  inequality  constraints. 
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